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

$C\frac{dV}{dt}=I.$

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:

(1)$\frac{dV}{dt}=\frac{-1}{\tau}V+\frac{1}{C}I.$

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:

$\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,

$\frac{df}{dt}=af+g$

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:

$\frac{d}{dt}y=Ay$

where

$\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$

and

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

$\frac{dV}{dt}=\frac{-1}{\tau}V+\frac{1}{C}I$

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

$\frac{dV}{dt}=\frac{-1}{Tau}V+\frac{1}{C}\iota$

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

$y_{t+h}=y(t+h)=e^{A(t+h)}y(0)=e^{Ah}e^{At}y(0)=e^{Ah}y(t)=e^{Ah}y_t.$

The complete update for the neuron can be written as

$y_{t+h}=e^{Ah}y_t + x_{t+h}$

where

$\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}$$ 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.