Integrating neural models using exact integration

The simple integrate-and fire model

For the simple integrate-and-fire model the voltage \(V\) is given as a solution of the equation:


This is just the 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_{\text{th}}\), at which point a delta function spike occurs.

A shortcoming of the simple integrate-and-fire model is that it implements no time-dependent memory. If the model receives a below-threshold signal at some time, it will retain that voltage boost until it fires again. This characteristic is not in line with observed neuronal behavior.

The leaky integrate-and fire model

In the leaky integrate-and-fire model, the memory problem is solved by adding a “leak” term \(\frac{-1}{R}V\) (\(R\) is the resistance and \(\tau=RC\)) to the membrane potential:

\[\frac{dV}{dt}=\frac{-1}{\tau}V+\frac{1}{C}I. :label: membrane\]

This reflects the diffusion of ions that occurs through the membrane when some equilibrium is not reached in the cell.

Solving a homogeneous linear differential equation

To solve membrane we start by looking at a simpler differential equation:

\[\frac{df}{dt}=af\text{, where } f:\mathbb{R}\to\mathbb{R} \text{ and } a\in\mathbb{R}.\]

Here the solution is given by \(f(t)=e^{at}\).

Solving a non-homogeneous linear differential equation

When you add another function \(g\) to the right hand side of our linear differential equation,


this is now a non-homogeneous differential equation. Things (can) become more complicated.

Solving it with variation of constants

This kind of differential equation is usually solved with “variation of constants” which gives us the following solution:

\[f(t)=e^{ct}\int_{0}^t g(s)e^{-cs}ds.\]

This is obviously not a particularly handy solution since calculating the integral in every step is very costly.

Solving it with exact integration

With exact integration, these costly computations can be avoided.

Restrictions to \(g\)

But only for certain functions \(g\)! I.e. if \(g\) satisfies (is a solution of):

\[\left(\frac{d}{dt}\right)^n g= \sum_{i=1}^{n}a_i\left(\frac{d}{dt}\right)^{i-1} g\]

for some \(n\in \mathbb{N}\) and a sequence \((a_i)_{i\in\mathbb{N}}\subset \mathbb{R}\).

For example this would be the case for \(g=\frac{e}{\tau_{syn}}t e^{-t/\tau_{\text{syn}}}\) (an alpha funciton), where \(\tau_{\text{syn}}\) is the rise time.

Reformulating the problem

The non-homogeneous differential equation is reformulated as a multidimensional homogeneous linear differential equation:



\[\begin{split}A=\begin{pmatrix} a_{n} & a_{n-1} & \cdots & \cdots & a_1 & 0 \\ 1 & 0 & \cdots & 0 & 0 & 0 \\ 0 & \ddots & \ddots & \vdots & \vdots & \vdots \\ \vdots & \ddots & \ddots & 0 & 0 & 0 \\ 0 & 0 & \ddots & 1 & 0 & 0 \\ 0 & 0 & \cdots & 0 & \frac{1}{C} & -\frac{1}{\tau} \\ \end{pmatrix}\end{split}\]

by choosing \(y_1,...,y_n\) canonically as:

\[\begin{split}\begin{align*} y_1 &= \left(\frac{d}{dt}\right)^{n-1}g\\ \vdots &= \vdots\\ y_{n-1} &= \frac{d}{dt}g\\ y_{n} &= g\\ y_{n+1} &= f. \end{align*}\end{split}\]

This makes ist very easy to determine the solution as

\[y(t)= e^{At}y_0\]


\[y_{t+h}=y(t+h)=e^{A(t+h)}\cdot y_0=e^{Ah}\cdot e^{At}\cdot y_0=e^{Ah}\cdot y_t.\]

This means that once we have calculated \(A\), propagation consists of multiplications only.

Example: The leaky integrate and fire model with alpha-function shaped inputs (iaf_psc_alpha)

The dynamics of the membrane potential \(V\) is given by:


where \(\tau\) is the membrane time constant and \(C\) is the capacitance. \(I\) is the sum of the synaptic currents and any external input:

Postsynaptic currents are alpha-shaped, i.e. the time course of the synaptic current \(\iota\) due to one incoming spike is

\[\iota (t)= \frac{e}{\tau_{syn}}t e^{-t/\tau_{\text{syn}}}.\]

The total input \(I\) to the neuron at a certain time \(t\) is the sum of all incoming spikes at all grid points in time \(t_i\le t\) plus an additional piecewise constant external input \(I_{\text{ext}}\):

\[I(t)=\sum_{i\in\mathbb{N}, t_i\le t }\sum_{k\in S_{t_i}}\hat{\iota}_k \frac{e}{\tau_{\text{syn}}}(t-t_i) e^{-(t-t_i)/\tau_{\text{syn}}}+I_{\text{ext}}\]

\(S_t\) is the set of indices that deliver a spike to the neuron at time \(t\), \(\tau_{\text{syn}}\) is the rise time and \(\iota_k\) represents the “weight” of synapse \(k\).

Exact integration for the iaf_psc_alpha model

First we make the substitutions:

\[\begin{split}\begin{align*} y_1 &= \frac{d}{dt}\iota+\frac{1}{\tau_{syn}}\iota \\ y_2 &= \iota \\ y_3 &= V \end{align*}\end{split}\]

for the equation


we get the homogeneous differential equation (for \(y=(y_1,y_2,y_3)^t\))

\[\begin{split}\frac{d}{dt}y= Ay= \begin{pmatrix} \frac{1}{\tau_{syn}}& 0 & 0\\ 1 & \frac{1}{\tau_{syn}} & 0\\ 0 & \frac{1}{C} & -\frac {1}{\tau} \end{pmatrix} y.\end{split}\]

The solution of this differential equation is given by \(y(t)=e^{At}y(0)\) and can be solved stepwise for a fixed time step \(h\):


The complete update for the neuron can be written as

\[y_{t+h}=e^{Ah}y_t + x_{t+h}\]


\[\begin{split}x_{t+h}+\begin{pmatrix}\frac{e}{\tau_{\text{syn}}}\\0\\0\end{pmatrix}\sum_{k\in S_{t+h}}\hat{\iota}_k\end{split}\]

as the linearity of the system permits the initial conditions for all spikes arriving at a given grid point to be lumped together in the term \(x_{t+h}\). \(S_{t+h}\) is the set of indices \(k\in 1,....,K\) of synapses that deliver a spike to the neuron at time \(t+h\).

The matrix \(e^{Ah}\) is constructed in the C++ implementation of the iaf_psc_alpha model in NEST.

Every matrix entry is calculated twice. For inhibitory postsynaptic inputs (with a time constant \(\tau_{syn_{in}}\)) and excitatory postsynaptic inputs (with a time constant \(\tau_{syn_{ex}}\)).

The update is performed here. The first multiplication evolves the external input. The others are the multiplication of the matrix \(e^{Ah}\) with \(y\) (for inhibitory and excitatory inputs).

If synaptic and membrane time constants become very close, \(\tau_m\approx \tau_{syn}\), the matrix \(e^{Ah}\) becomes numerically unstable. NEST handles this gracefully as described in the IAF Integration Singularity notebook.

For more information see [1].