The optimize package supports parameter estimation of stochastic generalized integrate and fire (gLIF) neurons from patch clamp data. In this document, we shall describe the following:
- Model equations.
- Object oriented implementation.
- Optimization routines.
The neurons estimated by the optimize packages consist of two components:
The following difference equation defines the subthreshold model:
See fit_neuron.optimize.threshold.StochasticThresh for the implemenation of this difference equation.
Whenever a spike occurs (as determined by the threshold equations), the state of the neuron is updated as follows:
This update rule is implemented for general cases of spike induced currents by fit_neuron.optimize.subthreshold.Voltage.spike().
The Optimize package allows two variations of spike triggered currents.
In [MS2011] the spike triggered currents take the form of a sum of step functions. This sum can be written as:
where the indicator functions are defined as follows:
By looking at the difference equation for the voltage, one sees that the parameter vector defines the shape of the spike triggered currents, and the parameters define the intervals during which the shape is constant.
In this case, the update rule simply appends the current spike to the spiking history :
The Gerstner adaptation currents are implemented in fit_neuron.optimize.sic_lib.StepSic.
The equations are written here in a form that matches the implemenation. The equations are written differently in [MS2011] but are perfectly equivalent up to a linear transformation of the parameter vector .
An alternative form of spike triggered currents is used in [MN2009] and consists of exponentially decaying currents with an additive reset. The equations are as follows:
and the reset equation, applied whenever the neuron spikes, is:
The Mihalas-Niebur adaptation currents are implemented in fit_neuron.optimize.sic_lib.ExpDecay_sic.
When calling the optimization routine fit_neuron.optimize.fit_gLIF.fit_neuron(), the user has the ability to specify any spike induced object he/she would like, as long as the user defines the class of the spike triggered current to inherit from the following abstract class: fit_neuron.optimize.sic_lib.SicBase.
The stochastic neuron has the following hazard rate:
where the parameters are the indicator variables (see above), and the parameters are probability currents which shall be referred to as voltage chasing currents. These currents give the stochastic spike emission process a component that adapts to the history of the voltage. The equations used for the voltage chasing currents are:
When the neuron spikes, the voltage chasing currents are set to the reset potential:
The hazard rate is computed at each time step and compared to a uniformly distributed random number to determine whether the neuron spikes here. This computation is performed by fit_neuron.optimize.threshold.StochasticThresh.update_X_arr().
The parameter estimation algorithm provided by the fit_neuron.optimize.fit_gLIF() function proceeds as follows:
The subthreshold parameters are obtained via linear regression to the observed voltage differences.
The equation that is solved is the following:
The value of that minimizes this expression is the parameter vector chosen for the subthreshold object fit_neuron.optimize.subthreshold.Voltage.
The values of used above represent the values in the recorded voltage traces. The values of the spike induced currents are computed based on the recorded voltage values and the recorded spike times. Hence the estimation process resembles maximum likelihood.
If no voltage nonlinearity is provided, or if it is set to None, the parameter vector will still correspond to the vector above but with the voltage nonlinearity skipped.
as computed by fit_neuron.optimize.threshold.StochasticThresh.update_X_arr().
The probability of there being a spike in a time increment is
The probability of observing spikes at indices in the spiking set , and no spikes at the indices in the complement of this set , is:
Taking logs of both sides, we obtain
which may be approximated up to a constant as:
The ‘th elements of the gradient of this function w.r.t. are:
The elements of the Hessian matrix are:
It is trivially seen that the Hessian matrix is negative definite. Hence the negative log likelihood is a convex function of the parameters and convex optimization techniques are applicable here. We use a Newton algorithm to update the values of the parameters:
The most computationally expensive step in this process is the computation of the gradients and hessians, which must be done at every step. Significant speedups can be achieved by distributing the computations of the gradients and the hessians to multiple processors. This is done in fit_neuron.optimize.threshold.par_calc_log_like_update().
|[RB2005]||Brette, Romain, and Wulfram Gerstner. “Adaptive exponential integrate-and-fire model as an effective description of neuronal activity.” Journal of neurophysiology 94.5 (2005): 3637-3642.|
|[MN2009]||Mihalas, Stefan, and Ernst Niebur. “A generalized linear integrate-and-fire neural model produces diverse spiking behaviors.” Neural computation 21.3 (2009): 704-718.|
|[MS2011]||(1, 2, 3) Mensi, Skander, et al. “Parameter extraction and classification of three cortical neuron types reveals two distinct adaptation mechanisms.” Journal of neurophysiology 107.6 (2012): 1756-1775.|