Contents

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:

- Subthreshold dynamics, which capture the dynamics of the membrane potential, its
dependence on the input current, and its dependence on spike times. The subthreshold
dynamics are implemented by
`fit_neuron.optimize.subthreshold.Voltage`. - Threshold dynamics, which predicts the probability of a spike occuring
depending on the current voltage and the past history of the neuron.
The threshold dynamics are implemented by
`fit_neuron.optimize.threshold.StochasticThresh`.

The following difference equation defines the subthreshold model:

where

- is the current time
- is the time step, typically corresponding to the time step in the recorded patch clamp data
- is the set of previous spike times
- The are the
*spike triggered currents*or*spike induced currents*which depend only on the spiking history of the neuron, and optionally on the current value of the voltage . - The function is a voltage nonlinearity function that allows the model to have an upward spike initiation. Commonly used voltage nonlinearities are the quadratic and exponential functions.

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:

where

- is the reset voltage.
- is the refractory period of the neuron.
- is an optional
*update rule*.

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`.

Note

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:

- Extract spikes and spike shapes from the raw data.
- Take the voltage traces with the spike shapes removed, and estimate the subthreshold parameters by linear regression of the derivative of the voltage.
- Simulate the model neuron using the same inputs as the raw data, and force the spikes to happen at the times the biological neuron was observed to spike. This will produce simulated voltage traces.
- Fit the threshold parameters such that these parameters maximize the log likelihood of the obseved spike trains being emitted by the simulated voltage traces.

The subthreshold parameters are obtained via linear regression to the observed voltage differences.

The equation that is solved is the following:

where

and

The value of that minimizes this expression is the parameter
vector chosen for the subthreshold object `fit_neuron.optimize.subthreshold.Voltage`.

Note

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.

Note

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.

The threshold parameters are obtained via max likelihood of the observed spike train. Following
along the lines of [MS2011] we may re-write the threshold equation in *Threshold Equations* as follows:

where

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. |