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 derivate 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.
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.
To solve (1) we start by looking at a simpler differential equation:
$\frac{df}{dt}=af$, where $f:\mathbb{R}\to\mathbb{R}$ and $a\in\mathbb{R}$
Here the solution is given by $f(t)=e^{at}$.
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.
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.
With exact integration these costly computations can be avoided.
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.
The non-homogeneous differential equation is reformulated as a multidimensional homogeneous linear differential equation:
$\frac{d}{dt}y=Ay$ where
$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}$
by choosing $y_1,...,y_n$ canonically as:
$\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*}$
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.
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$.
First we make the substitutions:
$y_1=\frac{d}{dt}\iota+\frac{1}{\tau_{syn}}\iota$
$y_2=\iota$
$y_3=V$
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$)
$\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$.
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
$x_{t+h}+\begin{pmatrix}\frac{e}{\tau_{\text{syn}}}\\0\\0\end{pmatrix}\sum_{k\in S_{t+h}}\hat{\iota}_k$
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 post synaptic inputs (with a time constant $\tau_{syn_{in}}$) and excitatory post synaptic inputs (with a time constant $\tau_{syn_{ex}}$).
And 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)
[1] RotterV S & Diesmann M (1999) Exact simulation of time-invariant linear systems with applications to neuronal modeling. Biologial Cybernetics 81:381-402; which you will find here
In [ ]: