1Integrating neural models using exact integration 2================================================= 3 4The simple integrate-and fire model 5~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 6 7For the simple integrate-and-fire model the voltage :math:`V` is given as a solution of the equation: 8 9.. math:: 10 C\frac{dV}{dt}=I. 11 12This is just the derivate of the law of capacitance :math:`Q=CV`. When an input current is applied, the membrane voltage increases with time until it reaches a constant threshold :math:`V_{\text{th}}`, at which point a delta function spike occurs. 13 14A 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. 15 16The leaky integrate-and fire model 17~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 18 19In the leaky integrate-and-fire model, the memory problem is solved by adding a "leak" term :math:`\frac{-1}{R}V` (:math:`R` is the resistance and :math:`\tau=RC`) to the membrane potential: 20 21.. math:: 22 \frac{dV}{dt}=\frac{-1}{\tau}V+\frac{1}{C}I. 23 :label: membrane 24 25This reflects the diffusion of ions that occurs through the membrane when some equilibrium is not reached in the cell. 26 27Solving a homogeneous linear differential equation 28~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 29 30To solve :math:numref:`membrane` we start by looking at a simpler differential equation: 31 32.. math:: 33 \frac{df}{dt}=af\text{, where } f:\mathbb{R}\to\mathbb{R} \text{ and } a\in\mathbb{R}. 34 35Here the solution is given by :math:`f(t)=e^{at}`. 36 37Solving a non-homogeneous linear differential equation 38~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 39When you add another function :math:`g` to the right hand side of our linear differential equation, 40 41.. math:: 42 \frac{df}{dt}=af+g 43 44this is now a non-homogeneous differential equation. Things (can) become more complicated. 45 46Solving it with variation of constants 47^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 48 49This kind of differential equation is usually solved with "variation of constants" which gives us the following solution: 50 51.. math:: 52 f(t)=e^{ct}\int_{0}^t g(s)e^{-cs}ds. 53 54This is obviously not a particularly handy solution since calculating the integral in every step is very costly. 55 56Solving it with exact integration 57^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 58 59With exact integration, these costly computations can be avoided. 60 61Restrictions to :math:`g` 62------------------------ 63But only for certain functions :math:`g`! I.e. if :math:`g` satisfies (is a solution of): 64 65.. math:: 66 \left(\frac{d}{dt}\right)^n g= \sum_{i=1}^{n}a_i\left(\frac{d}{dt}\right)^{i-1} g 67 68for some :math:`n\in \mathbb{N}` and a sequence :math:`(a_i)_{i\in\mathbb{N}}\subset \mathbb{R}`. 69 70For example this would be the case for :math:`g=\frac{e}{\tau_{syn}}t e^{-t/\tau_{\text{syn}}}` (an alpha funciton), where :math:`\tau_{\text{syn}}` is the rise time. 71 72Reformulating the problem 73^^^^^^^^^^^^^^^^^^^^^^^^^ 74 75The non-homogeneous differential equation is reformulated as a multidimensional homogeneous linear differential equation: 76 77.. math:: 78 \frac{d}{dt}y=Ay 79 80where 81 82.. math:: 83 A=\begin{pmatrix} 84 a_{n} & a_{n-1} & \cdots & \cdots & a_1 & 0 \\ 85 1 & 0 & \cdots & 0 & 0 & 0 \\ 86 0 & \ddots & \ddots & \vdots & \vdots & \vdots \\ 87 \vdots & \ddots & \ddots & 0 & 0 & 0 \\ 88 0 & 0 & \ddots & 1 & 0 & 0 \\ 89 0 & 0 & \cdots & 0 & \frac{1}{C} & -\frac{1}{\tau} \\ 90 \end{pmatrix} 91 92by choosing :math:`y_1,...,y_n` canonically as: 93 94.. math:: 95 \begin{align*} 96 y_1 &= \left(\frac{d}{dt}\right)^{n-1}g\\ 97 \vdots &= \vdots\\ 98 y_{n-1} &= \frac{d}{dt}g\\ 99 y_{n} &= g\\ 100 y_{n+1} &= f. 101 \end{align*} 102 103This makes ist very easy to determine the solution as 104 105.. math:: 106 y(t)= e^{At}y_0 107 108and 109 110.. math:: 111 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. 112 113This means that once we have calculated :math:`A`, propagation consists of multiplications only. 114 115Example: The leaky integrate and fire model with alpha-function shaped inputs (iaf_psc_alpha) 116^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 117 118The dynamics of the membrane potential :math:`V` is given by: 119 120.. math:: 121 \frac{dV}{dt}=\frac{-1}{\tau}V+\frac{1}{C}I 122 123where :math:`\tau` is the membrane time constant and :math:`C` is the capacitance. :math:`I` is the sum of the synaptic currents and any external input: 124 125Postsynaptic currents are alpha-shaped, i.e. the time course of the synaptic current :math:`\iota` due to one incoming spike is 126 127.. math:: 128 \iota (t)= \frac{e}{\tau_{syn}}t e^{-t/\tau_{\text{syn}}}. 129 130The total input :math:`I` to the neuron at a certain time :math:`t` is the sum of all incoming spikes at all grid points in time :math:`t_i\le t` plus an additional piecewise constant external input :math:`I_{\text{ext}}`: 131 132.. math:: 133 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}} 134 135:math:`S_t` is the set of indices that deliver a spike to the neuron at time :math:`t`, :math:`\tau_{\text{syn}}` is the rise time and :math:`\iota_k` represents the "weight" of synapse :math:`k`. 136 137Exact integration for the iaf_psc_alpha model 138--------------------------------------------- 139 140First we make the substitutions: 141 142.. math:: 143 \begin{align*} 144 y_1 &= \frac{d}{dt}\iota+\frac{1}{\tau_{syn}}\iota \\ 145 y_2 &= \iota \\ 146 y_3 &= V 147 \end{align*} 148 149for the equation 150 151.. math:: 152 \frac{dV}{dt}=\frac{-1}{Tau}V+\frac{1}{C}\iota 153 154we get the homogeneous differential equation (for :math:`y=(y_1,y_2,y_3)^t`) 155 156.. math:: 157 \frac{d}{dt}y= Ay= 158 \begin{pmatrix} 159 \frac{1}{\tau_{syn}}& 0 & 0\\ 160 1 & \frac{1}{\tau_{syn}} & 0\\ 161 0 & \frac{1}{C} & -\frac {1}{\tau} 162 \end{pmatrix} 163 y. 164 165The solution of this differential equation is given by :math:`y(t)=e^{At}y(0)` and can be solved stepwise for a fixed time step :math:`h`: 166 167.. math:: 168 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. 169 170The complete update for the neuron can be written as 171 172.. math:: 173 y_{t+h}=e^{Ah}y_t + x_{t+h} 174 175where 176 177.. math:: 178 x_{t+h}+\begin{pmatrix}\frac{e}{\tau_{\text{syn}}}\\0\\0\end{pmatrix}\sum_{k\in S_{t+h}}\hat{\iota}_k 179 180as 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 :math:`x_{t+h}`. :math:`S_{t+h}` is the set of indices :math:`k\in 1,....,K` of synapses that deliver a spike to the neuron at time :math:`t+h`. 181 182The matrix :math:`e^{Ah}` in the C++ implementation of the model in NEST is constructed `here <https://github.com/nest/nest-simulator/blob/b3fc263e073f46f0732c10efb34fcc90f3b6771c/models/iaf_psc_alpha.cpp#L243>`_. 183 184Every matrix entry is calculated twice. For inhibitory postsynaptic inputs (with a time constant :math:`\tau_{syn_{in}}`) and excitatory postsynaptic inputs (with a time constant :math:`\tau_{syn_{ex}}`). 185 186And the update is performed `here <https://github.com/nest/nest-simulator/blob/b3fc263e073f46f0732c10efb34fcc90f3b6771c/models/iaf_psc_alpha.cpp#L305>`_. The first multiplication evolves the external input. The others are the multiplication of the matrix :math:`e^{Ah}` with :math:`y`. (For inhibitory and excitatory inputs) 187 188References 189~~~~~~~~~~ 190 191.. [1] RotterV S & Diesmann M (1999) Exact simulation of time-invariant linear 192 systems with applications to neuronal modeling. Biologial Cybernetics 193 81:381-402. DOI: https://doi.org/10.1007/s004220050570 194