Neuron Models 3: Ensembles
Following on from my previous blog post on the exponential integrate and fire neuron model we will work to create a probabilistic view of the model. We present a numerical scheme to approximate this
This is the third blog post in a series - you can find the previous blog post here
Justifying Gaussian White Noise
We first begin with a small diversion, in the previous HH and EIF neuron firing examples we have assumed some sort of Gaussian white noise as an input signal. We briefly mentioned that this is a reasonable assumption but we will justify this in a bit more detail here. First we note that each neuron will typically take input signal from the order of 10,000 neurons. As such even in a low firing rate scheme a neuron will likely receive relatively large amount of input spikes. We can express this signal as:
$$I_p(t) = \sum^{N}_{i=0} J_{i} \sum_k \delta(t-t_i^k)$$
Where:
$N$ is the number of connected neurons
$J_{i}$ is the synaptic connection strength from neuron $i$
$t_i^k$ is the time of the kth spike recieved from neuron i
If we assume that these spikes arrive in an uncorrelated, memoryless fashion in the form of a Poisson process and that the connection strengths are suitably small: $\langle J_i \rangle \ll V_{Th} - V_{Re}$ (where angle brackets denote population average). Then we can apply a diffusion approximation:
$$I_p(t) = \sum^{N}_{i=0} J_{i} \sum_k \delta(t-t_i^k) \approx \mu + \sigma \xi(t)$$
Where:
$\mu = \langle J_i \rangle N \nu $
$\sigma = \langle J_i^2 \rangle N \nu $
$\nu$ is the mean firing rate over all connected neurons
$\xi(t)$ is a Gaussian white noise process
Of course as with all approximations this is subject to "small sample size" and $N$ needs to be suitably large.
Fokker-Planck
Recall that we specified the EIF model with Gaussian white noise as having dynamics: $$\tau \frac{dV_m}{dt} = (V_L - V_m) + \Delta_T e^{\left( \frac{V_m - V_T}{\Delta_T} \right)} + \sigma \sqrt{2 \tau}\xi_t $$
This is nothing more than an Ito process of the form: $$dX_t = \mu(X_t,t)dt + \sigma(X_t, t)dW_t $$ With standard Wiener process $W_t$. The Fokker-Planck equation gives us a probability distribution of this process $p(x,t)$ through the PDE: $$\frac{\partial}{\partial t} p(x,t) = -\frac{\partial}{\partial x}\left[\mu(x,t)p(x,t)\right] + \frac{\partial^2}{\partial x^2}\left[\frac{1}{2}\sigma^2(x,t)p(x,t)\right] $$ This formula can also be extended to higher dimensions in an obvious way. The derivation of this formula is fairly involved so not included in this blog post, most good textbooks on stochastic analysis should have a derivation for the interested reader.
In the case of the EIF model we can thus write down: $$ \frac{\partial p}{\partial t} = \frac{\sigma^2}{\tau}\frac{\partial^2p}{\partial V_m^2} + \frac{\partial}{\partial V_m} \left[ \frac{(V_m - V_L - \psi(V_m) )}{\tau} p(V_m,t) \right] $$ With $\psi(V_m)$ represnting the exponential firing term.
By the continuity equation we can write: $$ \frac{\partial p}{\partial t} = - \frac{\partial J}{\partial V_m} $$
Where $J$ represents the flux. By using this relation in the Fokker-Planck equation and integrating over voltage we get: $$ J(V_m, t) = - \frac{\sigma^2}{\tau}\frac{\partial p}{\partial V_m} - \frac{(V_m - V_L - \psi(V_m) )}{\tau} p(V_m,t) $$
We can also note that: $$J(V_{Re}^+,t) = J(V_{Re}^-, t) + r(t)$$
Where $V_{Re}^\pm$ represents the limit from above (+) or below (-) the reset voltage. the function $r(t)$ represents the average neuron firing rate. This is due to the implementaion of the voltage reset mechanism post spike. We can also note that for $V_m < V_{Re}$ we have $J(V_m, t) = 0$ and for $V_m > V_{Re}$ we have $J(V_m, t) = - r(t)$. We can then solve the flux equation to give: $$P(V_m, t) = \frac{r(t)\tau}{\sigma^2} \int_{max(V_m,V_{Re})}^{V_{Th}} exp \left( -\sigma^2 \int_{V_m}^u (x - V_L - \psi(x) )dx \right)du $$
Since the probability measure needs to integrate to 1, we can then write: $$r(t) = \left( \frac{\tau}{\sigma^2} \int_{-\infty}^{V_{Th}} \left( \int_{max(V_m,V_{Re})}^{V_{Th}} exp \left( -\sigma^2 \int_{V_m}^u (x - V_L - \psi(x) )dx \right)du \right) dV_m \right)^{-1} $$
(Note under the scheme presented there is no time dependence to any of these equations. Under time dependent signals we would have to be more careful and typically further approximations are made.)
So far we have not allowed for the refractory period, we have assumed that after reset the voltage trajectories continue as normal. Given we have chosen a deterministic refractory period we can just add this to the euqation above: $$r_{ref}(t) = \left( \frac{\tau}{\sigma^2} \int_{-\infty}^{V_{Th}} \left( \int_{max(V_m,V_{Re})}^{V_{Th}} exp \left( -\sigma^2 \int_{V_m}^u (x - V_L - \psi(x) )dx \right)du \right) dV_m + T_{Ref} \right)^{-1} $$
We can see that this integral will not give rise to an analytic solution in the case of EIF neurons. The forward Euler scheme we relied upon in the past will not perform well here. Instead we will use a slightly different numerical scheme.
Numerical Integration
(This is taken from Richardson [2007] - see references for further details)
Presented now is a numerical scheme for calculating the firing rate. Recall from above:
$$ J(V_m, t) = - r(t) \Theta(V - V_{Re}) = - \frac{\sigma^2}{\tau}\frac{\partial p}{\partial V_m} - \frac{(V_m - V_L - \psi(V_m) )}{\tau} p(V_m,t) $$
Where $\Theta(V)$ is the Heaviside step-function. Re-arranged this gives: $$- \frac{\partial p}{\partial V_m} = - \frac {\tau}{\sigma^2} r(t) \Theta(V - V_{Re}) + \sigma^{-2}(V_m - V_L - \psi(V_m)) p(V_m,t) $$
Which is of the form: $$ \frac{\partial p}{\partial V_m} = G(V_m)p(V_m) + H(V_m) $$
By applying a voltage discretization scheme: $V_k = V_{Lb} + k \Delta_V $ with $V_n = V_{Th}$ we can write down: $$ p(V_{k-1}) = p(V_k) e^{\int^{V_k}_{V_{k-1}} G(V)dV} + \int^{V_k}_{V_{k-1}} H(V) e^{\int^V_{V_{k-1}}G(U)dU} $$
We can approximate this as: $$ p(V_{k-1}) = p(V_k) e^{\Delta_V G(V_k)} + \Delta_V H(V_k) \left( \frac{e^{\Delta_V G(V_k)} - 1}{\Delta_V G(V_k)} \right) $$
Substituting back in the necessary formulae for $G$ and $H$ gives: $$ p(V_{k-1}) = p(V_k) e^{\Delta_V \sigma^{-2}(V_k - V_L - \psi(V_k)) } + \Delta_V \frac{\tau}{\sigma^2}r(t) \Theta(V_k - V_{Re})\left( \frac{e^{\Delta_V \sigma^{-2}(V_k - V_L - \psi(V_k))} - 1}{\Delta_V \sigma^{-2}(V_k - V_L - \psi(V_k))} \right) $$
However this still has unknown $r(t)$ in it. If we apply a transform: $q(V,t) = \frac{p(V,t)}{r(t)}$ then: $\sum q(V_k) = (r(t))^{-1}$ and: $$ q(V_{k-1}) = q(V_k) e^{\Delta_V \sigma^{-2}(V_k - V_L - \psi(V_k)) } + \Delta_V \frac{\tau}{\sigma^2} \Theta(V_k - V_{Re})\left( \frac{e^{\Delta_V \sigma^{-2}(V_k - V_L - \psi(V_k))} - 1}{\Delta_V \sigma^{-2}(V_k - V_L - \psi(V_k))} \right)$$
To simplify this expression we define functions $A$ and $B$ so that: $$ q(V_{k-1}) = q(V_k) A(V_k) + \Theta(V_k - V_{Re}) B(V_k) $$
And so we can calculate the firing rate. This scheme has a much better performance than an Euler scheme. We instantiate the scheme with $q(V_n) = 0$, we also select a value $V_{Lb}$ as a cut-off to stop iterating. An implementation of this method can be seen below:
# Evaluating the solution to the Fokker-Planck Equation to calculate the firing rate of an EIF neuron subject to Gaussian white noise
import numpy as np
# Set model parameters
# Membrane time constant tau (ms) and leak reversal potential VL (mV)
tau = 30
VL = -70
# Spike sharpness DelT (mV) and exponential potential threshold VT (mV)
DelT = 3
VT = -60
# Variation in gaussian noise sig
sig = 25
# Set voltage spike threshold Vth (mV), reset voltage Vr (mV) and refractory period Tref (ms)
Vth = 30
Vr = -70
Tref = 5
# Set up additional parameters for solving Fokker-Planck. DelV (mV) and VLb (mV)
DelV = 0.001
VLb = -100
Steps = int(np.ceil((Vth - VLb) / DelV))
q = np.zeros(Steps)
V = np.arange(Steps)*DelV + VLb
# For ease define function psi
def psi(V):
return DelT * np.exp((V - VT) / DelT)
def A(V):
return np.exp(DelV * sig**-2 *(V - VL - psi(V)))
def B(V):
if A(V) == 1.0:
return DelV * tau * sig**-2
else:
return DelV * tau * sig**-2 * (A(V) - 1) / np.log(A(V))
# Shut off numpy divide errors
np.seterr(divide='ignore')
for i in range(Steps -1, 0, -1):
if V[i] > Vr:
q[i-1] = q[i]*A(V[i]) + B(V[i])
else:
q[i-1] = q[i]*A(V[i])
r = 1/(q.sum()/1000000 + Tref/1000)
print("Firing rate:", round(r,1), "Hz")
We can modify the previous EIF firing code to estimate the firing rate, the results should be similar (note: for this I used 10m time steps, it is a slow running code!):
#collapse
# Implementation of a noisy EIF neuron using a forward Euler scheme
# Reduce N for quicker running code
import numpy as np
# Set seed for repeatability
np.random.seed(123)
# Set time step dt (ms) and number of steps N
dt = 0.001
N = 10000000
# Set model parameters
# Membrane time constant tau (ms) and leak reversal potential VL (mV)
tau = 30
VL = -70
# Spike sharpness DelT (mV) and exponential potential threshold VT (mV)
DelT = 3
VT = -60
# Variation in gaussian noise sig
sig = 25
# Set voltage spike threshold Vth (mV), reset voltage Vr (mV) and refractory period Tref (ms)
Vth = 30
Vr = -70
Tref = 5
# Set up voltage Vold (mV) and spike count Sp
Vold = Vr
Sp = 0
# Set up refractory period counter Tc (ms)
Tc = 0
for i in range(1, N):
if Tc > 0:
Vnew = Vr
Tc -= 1
else:
Vtemp = Vold + dt/tau*(VL - Vold) + DelT*dt/tau*np.exp((Vold - VT)/DelT) + sig*np.sqrt(2*dt/tau)*np.random.normal(0,1,1)
if Vtemp > Vth:
Vnew = Vr
Tc = np.ceil(Tref/dt)
Sp += 1
else:
Vnew = Vtemp
Vold = Vnew
print("Estimated firing rate:", round(Sp/(N*dt/1000),1), "Hz")
Which we can see is similar to the solution of the Fokker-Planck equation. In the limit $N \to \infty$ and decreasing the lattice sizes these approximations should become much closer.
Conclusion
We have seen that by using the Fokker-Planck framework we are able to calculate the mean firing rate of the EIF neuron. We can also notice that the numerical scheme to integrate the Fokker-Planck runs significantly faster than taking a Monte-Carlo approximation by simulating the EIF directly. We can also notice that the Fokker-Planck framework is easy to extend (e.g. to modulated noise or other applied signals) and further we can extend this to allow for connected networks of neurons (I may write an additional blog post on this in the future but will likely end up being quite similar to this one).
References
- https://neuronaldynamics.epfl.ch/online/Ch13.html - Online Neuronal Dynamics Textbook by Wulfram Gerstner, Werner M. Kistler, Richard Naud and Liam Paninski
- How Spike Generation Mechanisms Determine the Neuronal Response to Fluctuating Inputs - Nicolas Fourcaud-Trocme´, David Hansel, Carl van Vreeswijk, and Nicolas Brunel [2003]
- Firing-rate response of linear and nonlinear integrate-and-fire neurons to modulated current-based and conductance-based synaptic drive - Magnus J Richardson [2007]