Hodgkin–Huxley

The Hodgkin–Huxley model (H&H model)^{[8][9][10][11]} is a model of the relationship between the flow of ionic currents across the neuronal cell membrane and the membrane voltage of the cell.^{[8][9][10][11]} It consists of a set of nonlinear differential equations describing the behavior of ion channels that permeate the cell membrane of the squid giant axon. Hodgkin and Huxley were awarded the 1963 Nobel Prize in Physiology or Medicine for this work.

It is important to note the voltage-current relationship, with multiple voltage-dependent currents charging the cell membrane of capacity C_m

${\displaystyle C_{\mathrm {m} }{\frac {dV(t)}{dt}}=-\sum _{i}I_{i}(t,V).}$

The above equation is the time derivative of the law of capacitance, Q = CV where the change of the total charge must be explained as the sum over the currents. Each current is given by

${\displaystyle I(t,V)=g(t,V)\cdot (V-V_{\mathrm {eq} })}$

where g(t,V) is the conductance, or inverse resistance, which can be expanded in terms of its maximal conductance ḡ and the activation and inactivation fractions m and h, respectively, that determine how many ions can flow through available membrane channels. This expansion is given by

${\displaystyle g(t,V)={\bar {g}}\cdot m(t,V)^{p}\cdot h(t,V)^{q}}$

and our fractions follow the first-order kinetics

${\displaystyle {\frac {dm(t,V)}{dt}}={\frac {m_{\infty }(V)-m(t,V)}{\tau _{\mathrm {m} }(V)}}=\alpha _{\mathrm {m} }(V)\cdot (1-m)-\beta _{\mathrm {m} }(V)\cdot m}$

with similar dynamics for h, where we can use either τ and m_∞ or α and β to define our gate fractions.

The Hodgkin–Huxley model may be extended to include additional ionic currents. Typically, these include inward Ca²⁺ and Na⁺ input currents, as well as several varieties of K⁺ outward currents, including a "leak" current.

The result can be at the small end of 20 parameters which one must estimate or measure for an accurate model. In a model of a complex system of neurons, numerical integration of the equations are computationally expensive. Careful simplifications of the Hodgkin–Huxley model are therefore needed.

The model can be reduced to two dimensions thanks to the dynamic relations which can be established between the gating variables.^[13] it is also possible to extend it to take into account the evolution of the concentrations (considered fixed in the original model).^[14][15]

Perfect Integrate-and-fire

One of the earliest models of a neuron is the perfect integrate-and-fire model (also called non-leaky integrate-and-fire), first investigated in 1907 by Louis Lapicque.^[16] A neuron is represented by its membrane voltage V which evolves in time during stimulation with an input current I(t) according

${\displaystyle I(t)=C{\frac {dV(t)}{dt}}}$

which is just the time derivative of the law of capacitance, Q = CV. When an input current is applied, the membrane voltage increases with time until it reaches a constant threshold V_th, at which point a delta function spike occurs and the voltage is reset to its resting potential, after which the model continues to run. The firing frequency of the model thus increases linearly without bound as input current increases.

The model can be made more accurate by introducing a refractory period t_ref that limits the firing frequency of a neuron by preventing it from firing during that period. For constant input I(t)=I the threshold voltage is reached after an integration time t_int=CV_thr/I after start from zero. After a reset, the refractory period introduces a dead time so that the total time until the next firing is t_ref+t_int . The firing frequency is the inverse of the total inter-spike interval (including dead time). The firing frequency as a function of a constant input current is therefore

${\displaystyle \,\!f(I)={\frac {I}{C_{\mathrm {} }V_{\mathrm {th} }+t_{\mathrm {ref} }I}}.}$

A shortcoming of this model is that it describes neither adaptation nor leakage. If the model receives a below-threshold short current pulse at some time, it will retain that voltage boost forever - until another input later makes it fire. This characteristic is not in line with observed neuronal behavior. The following extensions make the integrate-and-fire model more plausible from a biological point of view.

Leaky integrate-and-fire

The leaky integrate-and-fire model which can be traced back to Louis Lapicque,^[16] contains, compared to the non-leaky integrate-and-fire model a "leak" term in the membrane potential equation, reflecting the diffusion of ions through the membrane. The model equation looks like^[1]

${\displaystyle C_{\mathrm {m} }{\frac {dV_{\mathrm {m} }(t)}{dt}}=I(t)-{\frac {V_{\mathrm {m} }(t)}{R_{\mathrm {m} }}}}$

A neuron is represented by an RC circuit with a threshold. Each input pulse (e.g. caused by a spike from a different neuron) causes a short current pulse. Voltage decays exponentially. If the threshold is reached an output spike is generated and the voltage is reset.

where V_m is the voltage across the cell membrane and R_m is the membrane resistance. (The non-leaky integrate-and-fire model is retrieved in the limit R_m to infinity, i.e. if the membrane is a perfect insulator). The model equation is valid for arbitrary time-dependent input until a threshold V_th is reached; thereafter the membrane potential is reset.

For constant input, the minimum input to reach the threshold is I_th = V_th / R_m. Assuming a reset to zero, the firing frequency thus looks like

${\displaystyle f(I)={\begin{cases}0,&I\leq I_{\mathrm {th} }\\\left[t_{\mathrm {ref} }-R_{\mathrm {m} }C_{\mathrm {m} }\log \left(1-{\tfrac {V_{\mathrm {th} }}{IR_{\mathrm {m} }}}\right)\right]^{-1},&I>I_{\mathrm {th} }\end{cases}}}$

which converges for large input currents to the previous leak-free model with the refractory period.^[17] The model can also be used for inhibitory neurons.^[18][19]

The biggest disadvantage of the Leaky integrate-and-fire neuron is that it does not contain neuronal adaptation so that it cannot describe an experimentally measured spike train in response to constant input current.^[20] This disadvantage is removed in generalized integrate-and-fire models that also contain one or several adaptation-variables and are able to predict spike times of cortical neurons under current injection to a high degree of accuracy.^[21][22][23]

Adaptive integrate-and-fire

Neuronal adaptation refers to the fact that even in the presence of a constant current injection into the soma, the intervals between output spikes increase. An adaptive integrate-and-fire neuron model combines the leaky integration of voltage V with one or several adaptation variables w_k (see Chapter 6.1. in the textbook Neuronal Dynamics^[27])

${\displaystyle \tau _{\mathrm {m} }{\frac {dV_{\mathrm {m} }(t)}{dt}}=RI(t)-[V_{\mathrm {m} }(t)-E_{\mathrm {m} }]-R\sum _{k}w_{k}}$

${\displaystyle \tau _{k}{\frac {dw_{k}(t)}{dt}}=-a_{k}[V_{\mathrm {m} }(t)-E_{\mathrm {m} }]-w_{k}+b_{k}\tau _{k}\sum _{f}\delta (t-t^{f})}$

where ${\displaystyle \tau _{m}}$ is the membrane time constant, w_k is the adaptation current number, with index k, ${\displaystyle \tau _{k}}$ is the time constant of adaptation current w_k, E_m is the resting potential and t^f is the firing time of the neuron and the Greek delta denotes the Dirac delta function. Whenever the voltage reaches the firing threshold the voltage is reset to a value V_r below the firing threshold. The reset value is one of the important parameters of the model. The simplest model of adaptation has only a single adaptation variable w and the sum over k is removed.^[28]

Spike times and subthreshold voltage of cortical neuron models can be predicted by generalized integrate-and-fire models such as the adaptive integrate-and-fire model, the adaptive exponential integrate-and-fire model, or the spike response model. In the example here, adaptation is implemented by a dynamic threshold which increases after each spike.^[22][23]

Integrate-and-fire neurons with one or several adaptation variables can account for a variety of neuronal firing patterns in response to constant stimulation, including adaptation, bursting, and initial bursting.^[24][25][26] Moreover, adaptive integrate-and-fire neurons with several adaptation variables are able to predict spike times of cortical neurons under time-dependent current injection into the soma.^[22][23]

Fractional-order leaky integrate-and-fire

Recent advances in computational and theoretical fractional calculus lead to a new form of model called Fractional-order leaky integrate-and-fire.^[29][30] An advantage of this model is that it can capture adaptation effects with a single variable. The model has the following form^[30]

${\displaystyle I(t)-{\frac {V_{\mathrm {m} }(t)}{R_{\mathrm {m} }}}=C_{\mathrm {m} }{\frac {d^{\alpha }V_{\mathrm {m} }(t)}{d^{\alpha }t}}}$

Once the voltage hits the threshold it is reset. Fractional integration has been used to account for neuronal adaptation in experimental data.^[29]

'Exponential integrate-and-fire' and 'adaptive exponential integrate-and-fire'

In the exponential integrate-and-fire model,^[33] spike generation is exponential, following the equation:

${\displaystyle {\frac {dV}{dt}}-{\frac {R}{\tau _{m}}}I(t)={\frac {1}{\tau _{m}}}\left[E_{m}-V+\Delta _{T}\exp \left({\frac {V-V_{T}}{\Delta _{T}}}\right)\right].}$

where ${\displaystyle V}$ is the membrane potential, ${\displaystyle V_{T}}$ is the intrinsic membrane potential threshold, ${\displaystyle \tau _{m}}$ is the membrane time constant, ${\displaystyle E_{m}}$is the resting potential, and ${\displaystyle \Delta _{T}}$ is the sharpness of action potential initiation, usually around 1 mV for cortical pyramidal neurons.^[31] Once the membrane potential crosses ${\displaystyle V_{T}}$, it diverges to infinity in finite time.^[34] In numerical simulation the integration is stopped if the membrane potential hits an arbitrary threshold (much larger than ${\displaystyle V_{T}}$) at which the membrane potential is reset to a value V_r . The voltage reset value V_r is one of the important parameters of the model. Importantly, the right-hand side of the above equation contains a nonlinearity that can be directly extracted from experimental data.^[31] In this sense the exponential nonlinearity is strongly supported by experimental evidence.

In the adaptive exponential integrate-and-fire neuron ^[32] the above exponential nonlinearity of the voltage equation is combined with an adaptation variable w

${\displaystyle \tau _{m}{\frac {dV}{dt}}=RI(t)+\left[E_{m}-V+\Delta _{T}\exp \left({\frac {V-V_{T}}{\Delta _{T}}}\right)\right]-Rw}$

${\displaystyle \tau {\frac {dw(t)}{dt}}=-a[V_{\mathrm {m} }(t)-E_{\mathrm {m} }]-w+b\tau \delta (t-t^{f})}$

Firing pattern of initial bursting in response to a step current input generated with the Adaptive exponential integrate-and-fire model. Other Firing patterns can also be generated.^[26]

where w denotes the adaptation current with time scale ${\displaystyle \tau }$. Important model parameters are the voltage reset value V_r, the intrinsic threshold ${\displaystyle V_{T}}$, the time constants ${\displaystyle \tau }$ and ${\displaystyle \tau _{m}}$ as well as the coupling parameters a and b. The adaptive exponential integrate-and-fire model inherits the experimentally derived voltage nonlinearity ^[31] of the exponential integrate-and-fire model. But going beyond this model, it can also account for a variety of neuronal firing patterns in response to constant stimulation, including adaptation, bursting, and initial bursting.^[26] However, since the adaptation is in the form of a current, aberrant hyperpolarization may appear. This problem was solved by expressing it as a conductance.^[35]

Noisy input model (diffusive noise)

A neuron embedded in a network receives spike input from other neurons. Since the spike arrival times are not controlled by an experimentalist they can be considered as stochastic. Thus a (potentially nonlinear) integrate-and-fire model with nonlinearity f(v) receives two inputs: an input ${\displaystyle I(t)}$ controlled by the experimentalists and a noisy input current ${\displaystyle I^{\rm {noise}}(t)}$ that describes the uncontrolled background input.

${\displaystyle \tau _{m}{\frac {dV}{dt}}=f(V)+RI(t)+RI^{\text{noise}}(t)}$

Stein's model^[39] is the special case of a leaky integrate-and-fire neuron and a stationary white noise current ${\displaystyle I^{\rm {noise}}(t)=\xi (t)}$ with mean zero and unit variance. In the subthreshold regime, these assumptions yield the equation of the Ornstein–Uhlenbeck process

${\displaystyle \tau _{m}{\frac {dV}{dt}}=[E_{m}-V]+RI(t)+R\xi (t)}$

However, in contrast to the standard Ornstein–Uhlenbeck process, the membrane voltage is reset whenever V hits the firing threshold V_th .^[39] Calculating the interval distribution of the Ornstein–Uhlenbeck model for constant input with threshold leads to a first-passage time problem.^[39][41] Stein's neuron model and variants thereof have been used to fit interspike interval distributions of spike trains from real neurons under constant input current.^[41]

In the mathematical literature, the above equation of the Ornstein–Uhlenbeck process is written in the form

${\displaystyle dV=[E_{m}-V+RI(t)]{\frac {dt}{\tau _{m}}}+\sigma \,dW}$

where ${\displaystyle \sigma }$ is the amplitude of the noise input and dW are increments of a Wiener process. For discrete-time implementations with time step dt the voltage updates are^[27]

${\displaystyle \Delta V=[E_{m}-V+RI(t)]{\frac {\Delta t}{\tau _{m}}}+\sigma {\sqrt {\tau _{m}}}y}$

where y is drawn from a Gaussian distribution with zero mean unit variance. The voltage is reset when it hits the firing threshold V_th .

The noisy input model can also be used in generalized integrate-and-fire models. For example, the exponential integrate-and-fire model with noisy input reads

${\displaystyle \tau _{m}{\frac {dV}{dt}}=E_{m}-V+\Delta _{T}\exp \left({\frac {V-V_{T}}{\Delta _{T}}}\right)+RI(t)+R\xi (t)}$

For constant deterministic input ${\displaystyle I(t)=I_{0}}$ it is possible to calculate the mean firing rate as a function of ${\displaystyle I_{0}}$.^[42] This is important because the frequency-current relation (f-I-curve) is often used by experimentalists to characterize a neuron. It is also the transfer function in

The leaky integrate-and-fire with noisy input has been widely used in the analysis of networks of spiking neurons.^[43] Noisy input is also called 'diffusive noise' because it leads to a diffusion of the subthreshold membrane potential around the noise-free trajectory (Johannesma,^[44] The theory of spiking neurons with noisy input is reviewed in Chapter 8.2 of the textbook Neuronal Dynamics.^[27]

Noisy output model (escape noise)

In deterministic integrate-and-fire models, a spike is generated if the membrane potential V(t) hits the threshold ${\displaystyle V_{th}}$. In noisy output models, the strict threshold is replaced by a noisy one as follows. At each moment in time t, a spike is generated stochastically with instantaneous stochastic intensity or 'escape rate' ^[27]

${\displaystyle \rho (t)=f(V(t)-V_{th})}$

that depends on the momentary difference between the membrane voltage V(t) and the threshold ${\displaystyle V_{th}}$.^[40] A common choice for the 'escape rate' ${\displaystyle f}$ (that is consistent with biological data^[22]) is

${\displaystyle f(V-V_{th})={\frac {1}{\tau _{0}}}\exp[\beta (V-V_{th}]}$

Stochastic spike generation (noisy output) depends on the momentary difference between the membrane potential V(t) and the threshold. The membrane potential V of the spike response model (SRM) has two contributions.^[45][46] First, input current I is filtered by a first filter k. Second the sequence of output spikes S(t) is filtered by a second filter η and fed back. The resulting membrane V(t) potential is used to generate output spikes by a stochastic process ρ(t) with an intensity that depends on the distance between membrane potential and threshold. The spike response model (SRM) is closely related to the Generalized Linear Model (GLM).^[47][48]

where ${\displaystyle \tau _{0}}$is a time constant that describes how quickly a spike is fired once the membrane potential reaches the threshold and ${\displaystyle \beta }$ is a sharpness parameter. For ${\displaystyle \beta \to \infty }$ the threshold becomes sharp and spike firing occurs deterministically at the moment when the membrane potential hits the threshold from below. The sharpness value found in experiments^[22] is ${\displaystyle 1/\beta \approx 4mV}$ which means that neuronal firing becomes non-negligible as soon the membrane potential is a few mV below the formal firing threshold.

The escape rate process via a soft threshold is reviewed in Chapter 9 of the textbook Neuronal Dynamics.^[27]

For models in discrete time, a spike is generated with probability

${\displaystyle P_{F}(t_{n})=F[V(t_{n})-V_{th}]}$

that depends on the momentary difference between the membrane voltage V at time ${\displaystyle t_{n}}$ and the threshold ${\displaystyle V_{th}}$.^[49] The function F is often taken as a standard sigmoidal ${\displaystyle F(x)=0.5[1+\tanh(\gamma x)]}$ with steepness parameter ${\displaystyle \gamma }$,^[40] similar to the update dynamics in artificial neural networks. But the functional form of F can also be derived from the stochastic intensity ${\displaystyle f}$ in continuous time introduced above as ${\displaystyle F(y_{n})\approx 1-\exp[y_{n}\Delta t]}$ where ${\displaystyle y_{n}=V(t_{n})-V_{th}}$ is the threshold distance.^[40]

Integrate-and-fire models with output noise can be used to predict the PSTH of real neurons under arbitrary time-dependent input.^[22] For non-adaptive integrate-and-fire neurons, the interval distribution under constant stimulation can be calculated from stationary renewal theory. ^[27]

Spike response model (SRM)

main article: Spike response model

The spike response model (SRM) is a general linear model for the subthreshold membrane voltage combined with a nonlinear output noise process for spike generation.^[40][52][50] The membrane voltage V(t) at time t is

${\displaystyle V(t)=\sum _{f}\eta (t-t^{f})+\int _{0}^{\infty }\kappa (s)I(t-s)\,ds+V_{\mathrm {rest} }}$

where t^f is the firing time of spike number f of the neuron, V_rest is the resting voltage in the absence of input, I(t-s) is the input current at time t-s and ${\displaystyle \kappa (s)}$ is a linear filter (also called kernel) that describes the contribution of an input current pulse at time t-s to the voltage at time t. The contributions to the voltage caused by a spike at time ${\displaystyle t^{f}}$ are described by the refractory kernel ${\displaystyle \eta (t-t^{f})}$. In particular, ${\displaystyle \eta (t-t^{f})}$ describes the reset after the spike and the time course of the spike-afterpotential following a spike. It therefore expresses the consequences of refractoriness and adaptation.^[40][23] The voltage V(t) can be interpreted as the result of an integration of the differential equation of a leaky integrate-and-fire model coupled to an arbitrary number of spike-triggered adaptation variables.^[24]

Spike firing is stochastic and happens with a time-dependent stochastic intensity (instantaneous rate)

${\displaystyle f(V-\vartheta (t))={\frac {1}{\tau _{0}}}\exp[\beta (V-\vartheta (t)]}$

with parameters ${\displaystyle \tau _{0}}$ and ${\displaystyle \beta }$ and a dynamic threshold ${\displaystyle \vartheta (t)}$ given by

${\displaystyle \vartheta (t)=\vartheta _{0}+\sum _{f}\theta _{1}(t-t^{f})}$

Here ${\displaystyle \vartheta _{0}}$ is the firing threshold of an inactive neuron and ${\displaystyle \theta _{1}(t-t^{f})}$ describes the increase of the threshold after a spike at time ${\displaystyle t^{f}}$.^[22][23] In case of a fixed threshold, one sets ${\displaystyle \theta _{1}(t-t^{f})}$=0. For ${\displaystyle \beta \to \infty }$ the threshold process is deterministic.^[27]

The time course of the filters ${\displaystyle \eta ,\kappa ,\theta _{1}}$ that characterize the spike response model can be directly extracted from experimental data.^[23] With optimized parameters the SRM describes the time course of the subthreshold membrane voltage for time-dependent input with a precision of 2mV and can predict the timing of most output spikes with a precision of 4ms.^[22][23] The SRM is closely related to linear-nonlinear-Poisson cascade models (also called Generalized Linear Model).^[48] The estimation of parameters of probabilistic neuron models such as the SRM using methods developed for Generalized Linear Models^[53] is discussed in Chapter 10 of the textbook Neuronal Dynamics.^[27]

Spike arrival causes postsynaptic potentials (red lines) which are summed. If the total voltage V reaches a threshold (dashed blue line) a spike is initiated (green) which also includes a spike-afterpotential. The threshold increases after each spike. Postsynaptic potentials are the response to incoming spikes while the spike-afterpotential is the response to outgoing spikes.

The name spike response model arises because in a network, the input current for neuron i is generated by the spikes of other neurons so that in the case of a network the voltage equation becomes

${\displaystyle V_{i}(t)=\sum _{f}\eta _{i}(t-t_{i}^{f})+\sum _{j=1}^{N}w_{ij}\sum _{f'}\varepsilon _{ij}(t-t_{j}^{f'})+V_{\mathrm {rest} }}$

where ${\displaystyle t_{j}^{f'}}$ is the firing times of neuron j (i.e., its spike train), and ${\displaystyle \eta _{i}(t-t_{i}^{f})}$ describes the time course of the spike and the spike after-potential for neuron i, ${\displaystyle w_{ij}}$ and ${\displaystyle \varepsilon _{ij}(t-t_{j}^{f'})}$ describe the amplitude and time course of an excitatory or inhibitory postsynaptic potential (PSP) caused by the spike ${\displaystyle t_{j}^{f'}}$ of the presynaptic neuron j. The time course ${\displaystyle \varepsilon _{ij}(s)}$ of the PSP results from the convolution of the postsynaptic current ${\displaystyle I(t)}$ caused by the arrival of a presynaptic spike from neuron j with the membrane filter ${\displaystyle \kappa (s)}$.^[27]

SRM0

The SRM0^[50][54][55] is a stochastic neuron model related to time-dependent nonlinear renewal theory and a simplification of the Spike Response Model (SRM). The main difference to the voltage equation of the SRM introduced above is that in the term containing the refractory kernel ${\displaystyle \eta (s)}$ there is no summation sign over past spikes: only the most recent spike (denoted as the time ${\displaystyle {\hat {t}}}$) matters. Another difference is that the threshold is constant. The model SRM0 can be formulated in discrete or continuous time. For example, in continuous time, the single-neuron equation is

${\displaystyle V(t)=\eta (t-{\hat {t}})+\int _{0}^{\infty }\kappa (s)I(t-s)\,ds+V_{\mathrm {rest} }}$

and the network equations of the SRM0 are^[50]

${\displaystyle V_{i}(t\mid {\hat {t}}_{i})=\eta _{i}(t-{\hat {t}}_{i})+\sum _{j}w_{ij}\sum _{f}\varepsilon _{ij}(t-{\hat {t}}_{i},t-t^{f})+V_{\mathrm {rest} }}$

where ${\displaystyle {\hat {t}}_{i}}$ is the last firing time neuron i. Note that the time course of the postsynaptic potential ${\displaystyle \varepsilon _{ij}}$ is also allowed to depend on the time since the last spike of neuron i to describe a change in membrane conductance during refractoriness.^[54] The instantaneous firing rate (stochastic intensity) is

${\displaystyle f(V-\vartheta )={\frac {1}{\tau _{0}}}\exp[\beta (V-V_{th})]}$

where ${\displaystyle V_{th}}$ is a fixed firing threshold. Thus spike firing of neuron i depends only on its input and the time since neuron i has fired its last spike.

With the SRM0, the interspike-interval distribution for constant input can be mathematically linked to the shape of the refractory kernel ${\displaystyle \eta }$ .^[40][50] Moreover the stationary frequency-current relation can be calculated from the escape rate in combination with the refractory kernel ${\displaystyle \eta }$.^[40][50] With an appropriate choice of the kernels, the SRM0 approximates the dynamics of the Hodgkin-Huxley model to a high degree of accuracy.^[54] Moreover, the PSTH response to arbitrary time-dependent input can be predicted.^[50]

Galves–Löcherbach model

3D visualization of the Galves–Löcherbach model for biological neural nets. This visualization is set for 4,000 neurons (4 layers with one population of inhibitory neurons and one population of excitatory neurons each) at 180 intervals of time.

The Galves–Löcherbach model^[56] is a stochastic neuron model closely related to the spike response model SRM0 ^[55][50] and the leaky integrate-and-fire model. It is inherently stochastic and, just like the SRM0 linked to time-dependent nonlinear renewal theory. Given the model specifications, the probability that a given neuron ${\displaystyle i}$ spikes in a period ${\displaystyle t}$ may be described by

${\displaystyle \mathop {\mathrm {Prob} } (X_{t}(i)=1\mid {\mathcal {F}}_{t-1})=\varphi _{i}{\Biggl (}\sum _{j\in I}W_{j\rightarrow i}\sum _{s=L_{t}^{i}}^{t-1}g_{j}(t-s)X_{s}(j),~~~t-L_{t}^{i}{\Biggl )},}$

where ${\displaystyle W_{j\rightarrow i}}$ is a synaptic weight, describing the influence of neuron ${\displaystyle j}$ on neuron ${\displaystyle i}$, ${\displaystyle g_{j}}$ expresses the leak, and ${\displaystyle L_{t}^{i}}$ provides the spiking history of neuron ${\displaystyle i}$ before ${\displaystyle t}$, according to

${\displaystyle L_{t}^{i}=\sup\{s<t:X_{s}(i)=1\}.}$

Importantly, the spike probability of neuron i depends only on its spike input (filtered with a kernel ${\displaystyle g_{j}}$ and weighted with a factor ${\displaystyle W_{j\to i}}$) and the timing of its most recent output spike (summarized by ${\displaystyle t-L_{t}^{i}}$).

FitzHugh–Nagumo

Sweeping simplifications to Hodgkin–Huxley were introduced by FitzHugh and Nagumo in 1961 and 1962. Seeking to describe "regenerative self-excitation" by a nonlinear positive-feedback membrane voltage and recovery by a linear negative-feedback gate voltage, they developed the model described by^[57]

${\displaystyle {\begin{array}{rcl}{\dfrac {dV}{dt}}&=&V-V^{3}/3-w+I_{\mathrm {ext} }\\\tau {\dfrac {dw}{dt}}&=&V-a-bw\end{array}}}$

where we again have a membrane-like voltage and input current with a slower general gate voltage w and experimentally-determined parameters a = -0.7, b = 0.8, τ = 1/0.08. Although not derivable from biology, the model allows for a simplified, immediately available dynamic, without being a trivial simplification.^[58] The experimental support is weak, but the model is useful as a didactic tool to introduce dynamics of spike generation through phase plane analysis. See Chapter 7 in the textbook Methods of Neuronal Modeling^[59]

Morris–Lecar

In 1981 Morris and Lecar combined the Hodgkin–Huxley and FitzHugh–Nagumo models into a voltage-gated calcium channel model with a delayed-rectifier potassium channel, represented by

${\displaystyle {\begin{aligned}&C{\frac {dV}{dt}}&=&-I_{\mathrm {ion} }(V,w)+I\\[6pt]&{\frac {dw}{dt}}&=&\varphi \cdot {\frac {w_{\infty }-w}{\tau _{w}}}\end{aligned}}}$

where ${\displaystyle I_{\mathrm {ion} }(V,w)={\bar {g}}_{\mathrm {Ca} }m_{\infty }\cdot (V-V_{\mathrm {Ca} })+{\bar {g}}_{\mathrm {K} }w\cdot (V-V_{\mathrm {K} })+{\bar {g}}_{\mathrm {L} }\cdot (V-V_{\mathrm {L} })}$.^[17] The experimental support of the model is weak, but the model is useful as a didactic tool to introduce dynamics of spike generation through phase plane analysis. See Chapter 7^[60] in the textbook Methods of Neuronal Modeling.^[59]

A two-dimensional neuron model very similar to the Morris-Lecar model can be derived step-by-step starting from the Hodgkin-Huxley model. See Chapter 4.2 in the textbook Neuronal Dynamics.^[27]

Hindmarsh–Rose

Building upon the FitzHugh–Nagumo model, Hindmarsh and Rose proposed in 1984^[61] a model of neuronal activity described by three coupled first-order differential equations:

${\displaystyle {\begin{aligned}&{\frac {dx}{dt}}&=&y+3x^{2}-x^{3}-z+I\\[6pt]&{\frac {dy}{dt}}&=&1-5x^{2}-y\\[6pt]&{\frac {dz}{dt}}&=&r\cdot (4(x+{\tfrac {8}{5}})-z)\end{aligned}}}$

with r² = x² + y² + z², and r ≈ 10⁻² so that the z variable only changes very slowly. This extra mathematical complexity allows a great variety of dynamic behaviors for the membrane potential, described by the x variable of the model, which includes chaotic dynamics. This makes the Hindmarsh–Rose neuron model very useful, because it is still simple, allows a good qualitative description of the many different firing patterns of the action potential, in particular bursting, observed in experiments. Nevertheless, it remains a toy model and has not been fitted to experimental data. It is widely used as a reference model for bursting dynamics.^[61]

Theta model and quadratic integrate-and-fire

The theta model, or Ermentrout–Kopell canonical Type I model, is mathematically equivalent to the quadratic integrate-and-fire model which in turn is an approximation to the exponential integrate-and-fire model and the Hodgkin-Huxley model. It is called a canonical model because it is one of the generic models for constant input close to the bifurcation point, which means close to the transition from silent to repetitive firing.^[62][63]

The standard formulation of the theta model is^[27][62][63]

${\displaystyle {\frac {d\theta (t)}{dt}}=(I-I_{0})[1+\cos(\theta )]+[1-\cos(\theta )]}$

The equation for the quadratic integrate-and-fire model is (see Chapter 5.3 in the textbook Neuronal Dynamics ^[27]))

${\displaystyle \tau _{\mathrm {m} }{\frac {dV_{\mathrm {m} }(t)}{dt}}=(I-I_{0})R+[V_{\mathrm {m} }(t)-E_{\mathrm {m} }][V_{\mathrm {m} }(t)-V_{\mathrm {T} }]}$

The equivalence of theta model and quadratic integrate-and-fire is for example reviewed in Chapter 4.1.2.2 of spiking neuron models.^[1]

For input I(t) that changes over time or is far away from the bifurcation point, it is preferable to work with the exponential integrate-and-fire model (if one wants the stay in the class of one-dimensional neuron models), because real neurons exhibit the nonlinearity of the exponential integrate-and-fire model.^[31]

The non-homogeneous Poisson process model (Siebert)

Siebert^[64][65] modeled the neuron spike firing pattern using a non-homogeneous Poisson process model, following experiments involving the auditory system.^[64][65] According to Siebert, the probability of a spiking event at the time interval ${\displaystyle [t,t+\Delta _{t}]}$ is proportional to a non-negative function ${\displaystyle g[s(t)]}$, where ${\displaystyle s(t)}$ is the raw stimulus.:

${\displaystyle P_{\text{spike}}(t\in [t',t'+\Delta _{t}])=\Delta _{t}\cdot g[s(t)]}$

Siebert considered several functions as ${\displaystyle g[s(t)]}$, including ${\displaystyle g[s(t)]\propto s^{2}(t)}$ for low stimulus intensities.

The main advantage of Siebert's model is its simplicity. The shortcomings of the model is its inability to reflect properly the following phenomena:

• The transient enhancement of the neuronal firing activity in response to a step stimulus.
• The saturation of the firing rate.
• The values of inter-spike-interval-histogram at short intervals values (close to zero).

These shortcomings are addressed by the age-dependent point process model and the two-state Markov Model.^[66][67][68]

Refractoriness and age-dependent point process model

Berry and Meister^[69] studied neuronal refractoriness using a stochastic model that predicts spikes as a product of two terms, a function f(s(t)) that depends on the time-dependent stimulus s(t) and one a recovery function ${\displaystyle w(t-{\hat {t}})}$ that depends on the time since the last spike

${\displaystyle \rho (t)=f(s(t))w(t-{\hat {t}})}$

The model is also called an inhomogeneous Markov interval (IMI) process.^[70] Similar models have been used for many years in auditory neuroscience.^[71][72][73] Since the model keeps memory of the last spike time it is non-Poisson and falls in the class of time-dependent renewal models.^[27] It is closely related to the model SRM0 with exponential escape rate.^[27] Importantly, it is possible to fit parameters of the age-dependent point process model so as to describe not just the PSTH response, but also the interspike-interval statistics.^[70][71][73]

Linear-nonlinear Poisson cascade model and GLM

The linear-nonlinear-Poisson cascade model is a cascade of a linear filtering process followed by a nonlinear spike generation step.^[74] In the case that output spikes feed back, via a linear filtering process, we arrive at a model that is known in the neurosciences as Generalized Linear Model (GLM).^[48][53] The GLM is mathematically equivalent to the spike response model SRM) with escape noise; but whereas in the SRM the internal variables are interpreted as the membrane potential and the firing threshold, in the GLM the internal variables are abstract quantities that summarizes the net effect of input (and recent output spikes) before spikes are generated in the final step.^[27][48]

The two-state Markov model (Nossenson & Messer)

The spiking neuron model by Nossenson & Messer^[66][67][68] produces the probability of the neuron firing a spike as a function of either an external or pharmacological stimulus.^[66][67][68] The model consists of a cascade of a receptor layer model and a spiking neuron model, as shown in Fig 4. The connection between the external stimulus to the spiking probability is made in two steps: First, a receptor cell model translates the raw external stimulus to neurotransmitter concentration, and then, a spiking neuron model connects neurotransmitter concentration to the firing rate (spiking probability). Thus, the spiking neuron model by itself depends on neurotransmitter concentration at the input stage.^[66][67][68]

Fig 4: High level block diagram of the receptor layer and neuron model by Nossenson & Messer.^[66][68]

Fig 5. The prediction for the firing rate in response to a pulse stimulus as given by the model by Nossenson & Messer.^[66][68]

An important feature of this model is the prediction for neurons firing rate pattern which captures, using a low number of free parameters, the characteristic edge emphasized response of neurons to a stimulus pulse, as shown in Fig. 5. The firing rate is identified both as a normalized probability for neural spike firing, and as a quantity proportional to the current of neurotransmitters released by the cell. The expression for the firing rate takes the following form:

${\displaystyle R_{\text{fire}}(t)={\frac {P_{\text{spike}}(t;\Delta _{t})}{\Delta _{t}}}=[y(t)+R_{0}]\cdot P_{0}(t)}$

where,

• P0 is the probability of the neuron being "armed" and ready to fire. It is given by the following differential equation:

${\displaystyle {\dot {P}}_{0}=-[y(t)+R_{0}+R_{1}]\cdot P_{0}(t)+R_{1}}$

P0 could be generally calculated recursively using the Euler method, but in the case of a pulse of stimulus, it yields a simple closed-form expression.^[66][75]

• y(t) is the input of the model and is interpreted as the neurotransmitter concentration on the cell surrounding (in most cases glutamate). For an external stimulus it can be estimated through the receptor layer model:

${\displaystyle y(t)\simeq g_{\text{gain}}\cdot \langle s^{2}(t)\rangle ,}$

with ${\displaystyle \langle s^{2}(t)\rangle }$ being a short temporal average of stimulus power (given in Watt or other energy per time unit).

• R₀ corresponds to the intrinsic spontaneous firing rate of the neuron.
• R₁ is the recovery rate of the neuron from the refractory state.

Other predictions by this model include:

1) The averaged evoked response potential (ERP) due to the population of many neurons in unfiltered measurements resembles the firing rate.^[68]

2) The voltage variance of activity due to multiple neuron activity resembles the firing rate (also known as Multi-Unit-Activity power or MUA).^[67][68]

3) The inter-spike-interval probability distribution takes the form a gamma-distribution like function.^[66][75]

Synaptic transmission (Koch & Segev)

According to the model by Koch and Segev,^[17] the response of a neuron to individual neurotransmitters can be modeled as an extension of the classical Hodgkin–Huxley model with both standard and nonstandard kinetic currents. Four neurotransmitters primarily influence the CNS. AMPA/kainate receptors are fast excitatory mediators while NMDA receptors mediate considerably slower currents. Fast inhibitory currents go through GABA_A receptors, while GABA_B receptors mediate by secondary G-protein-activated potassium channels. This range of mediation produces the following current dynamics:

• ${\displaystyle I_{\mathrm {AMPA} }(t,V)={\bar {g}}_{\mathrm {AMPA} }\cdot [O]\cdot (V(t)-E_{\mathrm {AMPA} })}$
• ${\displaystyle I_{\mathrm {NMDA} }(t,V)={\bar {g}}_{\mathrm {NMDA} }\cdot B(V)\cdot [O]\cdot (V(t)-E_{\mathrm {NMDA} })}$
• ${\displaystyle I_{\mathrm {GABA_{A}} }(t,V)={\bar {g}}_{\mathrm {GABA_{A}} }\cdot ([O_{1}]+[O_{2}])\cdot (V(t)-E_{\mathrm {Cl} })}$
• ${\displaystyle I_{\mathrm {GABA_{B}} }(t,V)={\bar {g}}_{\mathrm {GABA_{B}} }\cdot {\tfrac {[G]^{n}}{[G]^{n}+K_{\mathrm {d} }}}\cdot (V(t)-E_{\mathrm {K} })}$

where ḡ is the maximal^[8][17] conductance (around 1S) and E is the equilibrium potential of the given ion or transmitter (AMDA, NMDA, Cl, or K), while [O] describes the fraction of open receptors. For NMDA, there is a significant effect of magnesium block that depends sigmoidally on the concentration of intracellular magnesium by B(V). For GABA_B, [G] is the concentration of the G-protein, and K_d describes the dissociation of G in binding to the potassium gates.

The dynamics of this more complicated model have been well-studied experimentally and produce important results in terms of very quick synaptic potentiation and depression, that is fast, short-term learning.

The stochastic model by Nossenson and Messer translates neurotransmitter concentration at the input stage to the probability of releasing neurotransmitter at the output stage.^[66][67][68] For a more detailed description of this model, see the Two state Markov model section above.

The neurotransmitter-based energy detection scheme

The neurotransmitter-based energy detection scheme^[68][75] suggests that the neural tissue chemically executes a Radar-like detection procedure.

Fig. 6 The biological neural detection scheme as suggested by Nossenson et al.^[68][75]

As shown in Fig. 6, the key idea of the conjecture is to account for neurotransmitter concentration, neurotransmitter generation, and neurotransmitter removal rates as the important quantities in executing the detection task, while referring to the measured electrical potentials as a side effect that only in certain conditions coincide with the functional purpose of each step. The detection scheme is similar to a radar-like "energy detection" because it includes signal squaring, temporal summation, and a threshold switch mechanism, just like the energy detector, but it also includes a unit that emphasizes stimulus edges and a variable memory length (variable memory). According to this conjecture, the physiological equivalent of the energy test statistics is neurotransmitter concentration, and the firing rate corresponds to neurotransmitter current. The advantage of this interpretation is that it leads to a unit-consistent explanation which allows to bridge between electrophysiological measurements, biochemical measurements, and psychophysical results.

The evidence reviewed in^[68][75] suggest the following association between functionality to histological classification:

1. Stimulus squaring is likely to be performed by receptor cells.
2. Stimulus edge emphasizing and signal transduction is performed by neurons.
3. Temporal accumulation of neurotransmitters is performed by glial cells. Short-term neurotransmitter accumulation is likely to occur also in some types of neurons.
4. Logical switching is executed by glial cells, and it results from exceeding a threshold level of neurotransmitter concentration. This threshold crossing is also accompanied by a change in neurotransmitter leak rate.
5. Physical all-or-non movement switching is due to muscle cells and results from exceeding a certain neurotransmitter concentration threshold on muscle surroundings.

Note that although the electrophysiological signals in Fig.6 are often similar to the functional signal (signal power/neurotransmitter concentration / muscle force), there are some stages in which the electrical observation differs from the functional purpose of the corresponding step. In particular, Nossenson et al. suggested that glia threshold crossing has a completely different functional operation compared to the radiated electrophysiological signal and that the latter might only be a side effect of glia break.