The Leaky Integrate-and-Fire neuron model (LIF) was first described in 1907 by Louis Lapicque. This post derives the equations of the LIF neuron, compares exact and approximate methods to solve them, and simulates the Brunel network of excitatory and inhibitory LIF neurons.

The LIF model simplifies the complexity of neurons and assumes the neuron structure to be a single leaky compartment (point neuron model) with a membrane described by an RC circuit, as:

A neuron receives an input current \(I(t)\). Its membrane acts as a capacitor \(C\) in parallel with a resistor \(R\), in series with a battery of potential \(u_{\text{rest}}\). Adapted from: Neuronal Dynamics, Fig. 1.6

The main formulation of its electrical activity is given by the current equation, with the following terms:

  • The current in the capacitor \(C\) is given by \(I_C = \frac{dQ}{dt} = \frac{d CV}{dt} = C \frac{dV}{dt}\), for charge \(Q\) and capacitance \(C\);
  • The current in the resistance \(R\) is given by Ohm’s law as \(I_R = V/R\), for voltage \(V\) (relative to the resting potential \(u_{\text{rest}}\)) and resistance \(R\).

The law of conservation of electric charge states that the total charge of an isolated system remains constant. Thus, the currents entering a node of the circuit sum to the currents leaving it (Kirchhoff’s current law), i.e. the input current is \(I = I_C + I_R\). Equivalently:

\[C \frac{dV}{dt}(t) = -\frac{V(t)}{R} + I(t)\]

It is common to multiply all terms by \(R\) to expose the time constant \(\tau = R C\). This is a very important metric in the analysis of the behaviour of models, as it tells how fast the dynamics of the neuron are. The formulation then becomes:

\[\tau \frac{dV}{dt}(t) = -V(t) + RI(t) \label{equation_diff}\]

The external currents (such as synaptic activity) are added to \(I(t)\). In practice, the network contribution \(RI(t)\) is the sum of all currents delivered from any neuron \(j\) to a given neuron \(i\) over time, and is described as:

\[RI(t)_i = \tau \sum_{j} J_{ij} \sum_{t^\prime} \delta (t-t^\prime _{j} -D) \label{equation_current}\]

where \(J_{ij}\) is the postsynaptic potential (PSP) amplitude, \(\delta (t)\) the Dirac function, \(t^\prime _{j}\) the spiking times \(t^\prime\) of the pre-synaptic neuron \(j\), and \(D\) the transmission delay. For simplicity, \(V(t)\) is represented as \(V_t\).

The neuron’s potential is considered to rise steeply when there is an input current (by \(J_{ij}\) at each input spike), thus discarding the exponential charging of the capacitor. On the other hand, the discharge follows an exponential decay with the time constant \(\tau\).

To simulate the Action Potential (or spike) of biological neurons, the model assumes that each neuron has a certain fixed firing threshold. When the potential reaches the firing threshold, the neuron spikes (or fires), and its potential immediately drops to a reset potential. Axonal branches of neurons are simulated by a synaptic delay that regulates the time a spike signal takes to reach the connected post-synaptic neurons. This propagation of spike and current, which dictates the synaptic delay, is illustrated in the following picture:

An action potential travelling along an axon, in the direction of the arrows. The membrane is positive outside and negative inside, except at the action potential, where the polarity is reversed. Source: Action Potential.gif, Wikimedia Commons

Afterwards, for a given refractory period, the neuron potential stays constant at the reset value, to simulate the refraction of biological neurons. This may not be a very detailed simulation of neuron activity, but it is a good starting point for a simple neuron model, before moving on to more complex models in future posts.

Resolution

Two different stepping methods are derived from the analytical solution of the first-order ODE, both based on a two-step algorithm:

  • Calculate the decay from the previous step: \(V_t \leftarrow V_{t- \Delta t} \exp(\frac{-\Delta t}{\tau})\)
  • Add the voltage change due to the network current \(RI\) between \(t-\Delta t\) and \(t\): \(V_t \leftarrow V_t + RI_{(t-\Delta t, t]} / \tau\), where \(RI_{(t-\Delta t, t]}\) is the integral of \(RI\) over that interval.

Neurons can be advanced with fixed-step or variable-step integration:

  • Fixed step integration assumes the time step \(\Delta t\) to be constant throughout the execution.
  • Variable step methods advance the neuron with the smallest \(t\) in the network, with a \(\Delta t\) computed as the minimum of the following two values:
    • the time difference to the next incoming spike, so that each spike is added at its exact time;
    • \(t^\star + D - t\), where \(t^\star\) is the time of the neuron with the second smallest \(t\) and \(D\) is the synaptic delay, i.e. the largest step that can be taken from the neuron at \(t\) without missing a spike from the neuron at time \(t^\star\) (if it spikes).

The fixed and variable step integrations can also be computed numerically, instead of using the analytical solution. This is common practice for efficiency purposes, as the required exponential and division operations are computationally costly when \(\Delta t\) changes at every step. The solution is not exact, but it approximates the exact solution well enough when the timestep \(\Delta t\) is small enough. As a simple example, two approximate fixed-step Euler methods can be used:

  • Explicit Forward Euler:
\[\begin{align*} & \tau \frac{dV_{t - \Delta t}}{dt} = -V_{t - \Delta t} + RI_{(t-\Delta t, t]} / \Delta t \\ & V_t = V_{t - \Delta t} + \frac{dV_{t- \Delta t}}{dt} \Delta t \\ & \qquad \text{(forward difference)} \\ \Leftrightarrow & V_t = V_{t - \Delta t} + \frac{-V_{t - \Delta t} + RI_{(t-\Delta t, t]} / \Delta t}{\tau} \Delta t \end{align*}\]
  • Implicit Backward Euler:
\[\begin{align*} & \tau \frac{dV_{t}}{dt} = -V_{t} + RI_{(t-\Delta t, t]} / \Delta t \\ & \tau \frac{V_{t} - V_{t - \Delta t}}{\Delta t} = -V_{t} + RI_{(t-\Delta t, t]} / \Delta t \\ & \qquad \text{(backward difference)} \\ \Leftrightarrow & V_t = \frac{RI_{(t-\Delta t, t]} + \tau * V_{t-\Delta t}}{\Delta t + \tau} \end{align*}\]

with \(\frac{dV_{t}}{dt}\) and \(RI_{(t-\Delta t, t]}\) computed from equations \ref{equation_diff} and \ref{equation_current}, respectively. Note that \(RI_{(t-\Delta t, t]}\) is divided by \(\Delta t\) so that it is expressed per time unit instead of per timestep.

With a fixed step, however, the three methods have the same cost, as the factor \(\exp(\frac{-\Delta t}{\tau})\) of the analytical solution is a constant, computed only once. Explicit Euler is stable only for \(\Delta t < 2 \tau\), while implicit Euler is stable for any \(\Delta t\).

The variable step method is an event-driven algorithm: it is the most accurate, as each spike is added at its exact time, but it advances one neuron at a time. The fixed step methods are clock-driven: they advance all neurons together at every step, but align the spike times to the time grid (Brette et al., 2007). In practice, the fixed step suits large networks, as all neurons run the same operations at every step, while the variable step suits models where the exact timing of spikes matters.

For the sake of comparison, the picture below illustrates the potential over time of a given neuron in a network, for the aforementioned methods, on a \(100ms\) simulation of a (Brunel) network of \(10000\) neurons with default neuron parameters. Results are given for the Implicit and Explicit Euler methods with \(\Delta t = 0.1ms\), and for the analytical solution with fixed (\(\Delta t=0.1ms\)) and variable timestepping:

Note the vertical straight lines, which refer to the voltage reset after a spike, and the horizontal straight lines, which refer to the refractory period after firing. The four solutions diverge after the first spikes of the network, at around 14 ms, as a small difference in a spike time changes the input of all post-synaptic neurons.

Brunel Network

Network activity was mentioned in the previous section. Indeed, a neuron without any kind of stimulus (either electric current injection or synaptic currents) will always remain at rest.

The most common model of network activity is called the Brunel Network Model (Dynamics of Sparsely Connected Networks of Excitatory and Inhibitory Spiking Neurons, Brunel, Journal of Computational Neuroscience, 2000, pdf). It is highly relevant, as it provides a good estimation of network dynamics.

The model proposes a network of 10000 simulated LIF neurons (80% excitatory, 20% inhibitory; the paper simulates 12500), with external stimuli from excitatory neurons outside the network, which connect to both the excitatory and the inhibitory neurons. Excitatory means that the PSP is positive; the converse holds for inhibitory neurons. The external neurons are not simulated, so their spiking is generated by a Poisson random generator (PG). The diagram is the following:

The Brunel network: \(N_E\) excitatory neurons (E) and \(N_I\) inhibitory neurons (I). Each neuron receives \(C_E\) connections from excitatory neurons, \(C_I\) from inhibitory neurons, and \(C_E\) from Poisson generators (PG). Arrows are excitatory connections, and flat bars inhibitory ones. Adapted from: NEST by Example, Fig. 2

The tricky part here is to define the parameters of the network in terms of connectivity, synaptic weights, firing threshold, etc. The original paper describes them in detail, but for simplicity, they are summarized in the tables below:

Network Parameters
ParameterValue
Number of excitatory neurons \(N_E\)8000
Number of inhibitory neurons \(N_I\)2000
Excitatory synapses per neuron \(C_E\)800
Inhibitory synapses per neuron \(C_I\)200
Neuron Parameters
ParameterValue
Membrane time constant \(\tau\)20 ms
Refractory period \(\tau_{\text{rp}}\)2 ms
Firing threshold \(V_{\text{th}}\)20 mV
Membrane capacitance \(C\)1 pF
Resting potential \(u_{\text{rest}}\)0 mV
Reset potential \(V_{\text{reset}}\)10 mV
Excitatory PSP amplitude \(J_E\)0.1 mV
Inhibitory PSP amplitude \(J_I\)-0.5 mV
Synaptic delay \(D\)1.5 ms
Background rate \(\eta\)2.0

source: Marc-Oliver Gewaltig, ResearchGate

The stimuli parameters can now be adjusted to regulate network activity. The background rate \(\eta\) is the rate of the Poisson generators, in units of the rate that brings a neuron to its threshold without the input of the network. The relative strength \(g\) of the inhibitory synapses sets \(J_I = -g J_E\). Brunel defines a model for slow, asynchronous and fast network dynamics, which exhibit the following average neuron spiking rates (Hz) and global oscillation frequencies (Hz) of the network activity, where SI stands for synchronous irregular and AI for asynchronous irregular:

Firing rateGlobal frequency
SimulationTheorySimulationTheory
SI, fast60.7 Hz55.8 Hz180 Hz190 Hz
AI37.7 Hz38.0 Hz——
SI, slow5.5 Hz6.5 Hz22 Hz29 Hz

source: Brunel 2000

Notice how close the simulated network is to the theoretical predictions. For that reason, although it is a simple approximation, the Brunel Network is a powerful tool for studying network dynamics. The model is particularly useful in very long simulations of large networks, as it is computationally inexpensive compared to more detailed models such as the Hodgkin-Huxley model.

Output

A C++-based Brunel network model of LIF neurons was implemented, and the appropriate stimuli were injected to reproduce the results of the original paper. The outputs match the theoretical model in terms of spike rate and global oscillation frequency. Results are presented as a timeline (left to right) of neuron spike times (a.k.a. spike trace, top) and the number of neurons that spike at a given instant (bottom).

Slow dynamics (SI slow):

Asynchronous dynamics (AI):

Synchronous regular dynamics (SR):

In the pictures above, the global activity oscillates slowly in the SI slow state, where each neuron fires irregularly at a very low rate, and is stationary in the AI state. In the SR state, the network is almost fully synchronized, and the neurons fire regularly at high rates.

Other resources