# 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:

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 (1) we start by looking at a simpler differential equation:

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:

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):

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:

where

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

This makes ist very easy to determine the solution as

and

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

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}}\):

\(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:

for the equation

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

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

where

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}\) in the C++ implementation of the model in NEST is constructed here.

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.