Problem Set 3: Spike trains
AT2 – Neuromodeling: Problem set #3 SPIKE TRAINS
Younesse Kaddar
PROBLEM 1: Poisson spike trains
Link of the iPython notebook for the code
Brain neuron emit spikes seemingly randomly: we will aim to model such spike trains, by relying on the following assumptions:
 spike events occur independently of one another
 they occur at a constant firing rate $r$
 there exists a time bin $Δt$ within which no two spikes can be fired

within each time span $Δt$, the probability of the neuron firing is $rΔt$.

In other words: the random variable $X_t$ whose outcome is $1$ (resp. $0$) if there is a (resp. no) spike between $t$ and $t+Δt$ follows a Bernoulli distribution:
\[X_t \sim ℬ(r Δt)\]

1. Generating spike trains
Let $M ≝ T/Δt$, where $T$ is the length (in seconds) of the spike trains we will consider from now on.
In compliance with what we assumed before, a spike train can be modeled as a $M$dimensional random vector whose $i$th coordinate (for $0 ≤ i ≤ M1$) is the output of $X_{i Δt}$ (the values of which are in $\lbrace 0, 1 \rbrace$). In concrete terms, when implementing such a vector of $\lbrace 0, 1 \rbrace^M$: at each coordinate, it has a constant probability to be equal to $1$ ($0$ otherwise).
Example: For instance, if on average every fourth element in the vector is made to indicate a spike (i.e. is equal to $1$), that is: the previously mentioned constant probability is equal to $1/4$:
To clarify the time units, let us set
 $Δt ≝ 2 \text{ ms}$
 $T = 1 \text{ s}$
 $r = 25 \text{ spikes/sec}$
Then here is an example of a spike train:
2. Spike counts & Interspike intervals
Now, informative features can be extracted from each spike train, such as:
 its total number of spikes (which gives a global information about the spike train, but have us overlook the time step at which each spike occurred)
 the interspike intervals (which give us local clues on the contrary, but encompass the timerelated information we lost with the spike counts only)
A natural question to ask ourselves is: what is their distribution? To this end, let us generate $N ≝ 200$ spike trains with firing rate $r ≝ 25 \text{ Hz}$.
Here is the rastergram of $50$ of them:
Spike counts
Here is the spike count histogram of all of these $200$ spike trains:
We claim that spike counts are following a Poisson distribution of mean $rT ≝ 25$ (spikes), which may not seem obvious at first sight with so few spike trains:
However, it becomes more apparent when the number $N$ of spike trains is increased. For $N ≝ 100000$ for instance:
the spike counts distribution matches very well the Poisson one:
This can be explained theoretically: the assumptions stated in the introduction make the Poisson distribution perfectly suited to model the spike counts distribution. More precisely: by definition, the random variable $X$ counting the number of spikes is binomial:
\[X ≝ \sum\limits_{ n=0 }^{M1} \underbrace{X_{nΔt}}_{\; \sim \; ℬ(rΔt)} \; \sim \; ℬ(M, rΔt)\]as the $X_t$ are independent (by hypothesis) and identically distributed ($∀t, X_t \sim \; ℬ(rΔt)$).
But it is well known that the Poisson distribution $𝒫(\underbrace{M × rΔt}_{= \, rT})$ successfully approximates such a binomial distribution provided that $M$ is large enough and $rΔt$ is small enough. Here:
 $M ≝ T/Δt = \frac{1000}{2} = 500$
 $rΔt = 25 × 2 \cdot 10^{3} = 0.05$
which makes the Poisson distribution $𝒫(rT)$, i.e. $𝒫(25)$, perfectly suited to approximate the binomial distribution $ℬ(M, rΔt)$!
Interspike intervals
Here is the interspike interval histogram of $200$ spike trains with firing rate $r ≝ 25 \text{ Hz}$:
From a theoretical standpoint: if we denote by $Y$ the random variable whose outcome is the length of an interspike interval, in compliance with the assumptions in the introduction, it comes that:
\[\begin{cases} ∀t ≤ 0, \quad P(Y > t) = 1 \\ ∀t > 0, \quad P(Y > t) = (1  rΔt)^{t/Δt} = \exp\left(\frac{t}{Δt} \ln(1  rΔt)\right) = \exp\left( t \, \frac{\ln\left(\frac{1}{1  rΔt}\right)}{Δt}\right) \end{cases}\]So:
\[\begin{cases} ∀t ≤ 0, \quad P(Y ≤ t) = 0 \\ ∀t > 0, \quad P(Y ≤ t) = 1  \exp\left( t \, \ln\left(\frac{1}{1  rΔt}\right)/Δt\right) \end{cases}\]Therefore, by definition of its cumulative distribution function, the interspike interval variable $Y$ follows an exponential distribution of rate $\ln\left(\frac{1}{1  rΔt}\right)/Δt$
PROBLEM 2: Analysis of spike train
Link of the iPython notebook for the code
One records and analyses several spike trains from a single neuron in the primary somatosensory cortex of a monkey experiencing a vibratory stimulus on the fingertip, for different vibration frequencies $f_1 ∈ \lbrace 8.4, 12, 15.7, 19.6, 23.6, 25.9, 27.7, 35\rbrace \text{ Hz}$. The simulation is conducted between $200$ and $700$ msec.
As shown in Figure a. and Figure b., the spikes within the simulation period (between $200$ and $700$ msec):
 seem to be periodically spaced for a given frequency $f_1$
 are more and more numerous as the frequency increases
and even more: the number of spikes seem to depend linearly on the vibration frequency.
The latter observation is backed up by the shape of the average spike count, plotted as a function of the vibration frequency:
A linear regression of the average spike counts (as a function of the vibration frequency) can me made:
 Linear equation: \(y = 1.39 x + 3.08\)
 Correlation coefficient: $0.997$
 Standard error of the estimate: $0.04$
The very high correlation coefficient and very low standard error of the estimate highlight an almost perfect match!
We have similar results for the average firing rate as a function of the vibration frequency:
The linear regression yields:
 Linear equation: \(y = 2.78 x + 6.15\)
 Correlation coefficient: $0.997$
 Standard error of the estimate: $0.08$
As a result: the firing rate of the neuron can be regarded as a nondecreasing linear function of the vibration frequency.
PROBLEM 3: IntegrateandFire neuron
Link of the iPython notebook for the code
We will now:

first study the integrateandfire model with a constant input current, to acount for the creation of action potential (spikes) in real neurons

secondly make the input current vary over time, to try to account for the spike trains that were experimentally measured in the primary somatosensory cortex of a monkey experiencing a vibratory stimulus on the finger (cf. the previous problem)
1. IntegrateandFire neuron with constant input current
Voltage accross a passive membrane
The voltage across a passive neuron’s membrane when a current $I$ is injected is given by:
\[C \frac{ {\rm d}V}{ {\rm d}t} = g_L (E_L  V(t))+I\]where
 $C = 1 \text{ nF}$ is the membrane capacitance
 $g_L = 0.1 \text{ μS}$ is the leak conductance of the membrane
 and $E_L = −70 \text{ mV}$ its reversal potential
On top of that, we will set
 $I$ to be equal to $1 \text{ nA}$
 $V(0)$ to be equal to $E_L$
The Euler method yields:
\[\begin{align*} & \; C \frac{ ΔV}{ Δt} = g_L (E_L  V(t))+I\\ ⟺ & \; \frac{ ΔV}{ Δt} = \frac 1 C (g_L (E_L  V(t))+I)\\ ⟺ & \; V(t+Δt)  V(t) = \frac {Δt} C (g_L (E_L  V(t))+I)\\ ⟺ & \; V(t+Δt) = \Big(1 \frac {Δt \cdot g_L} C\Big) V(t) + \frac {Δt \, (g_L \, E_L + I)} C \end{align*}\]which enables us, by setting $∆t = 1 \text{ ms}$, to plot the voltage until time $t = 100 \text{ ms}$ for instance:
Now, for various current inputs $I$:
For positive currents: It appears that the voltage accross the passive membrane is a nondecreasing monotonous function of the time (reaching a plateau after $40$ ms in our case) which increases when the input current does so.
For the negative currents: the curves are symmetric to the positive current ones with respect to the $V = E_L$ line.
To be more precise, let us solve for $V$ analytically:
\[\begin{align*} & \; \frac{ {\rm d}V}{ {\rm d}t} = \frac {g_L} C (E_L  V(t))+I/C\\ ⟺ & \; \frac{ {\rm d}V}{ {\rm d}t} + \frac {g_L} C V(t) = \frac {g_L \, E_L +I} C\\ ⟺ & \; \frac{ {\rm d}V}{ {\rm d}t} \, {\rm e}^{\frac {g_L \, t} C} + \frac {g_L} C V(t) \, {\rm e}^{\frac {g_L \, t} C} = \frac {g_L \, E_L +I} C \, {\rm e}^{\frac {g_L \, t} C}\\ ⟺ & \; \frac{ {\rm d}}{ {\rm d}t} \left[V(t) \, {\rm e}^{\frac {g_L \, t} C}\right] = \frac {g_L \, E_L +I} C \, {\rm e}^{\frac {g_L \, t} C}\\ ⟺ & \; V(t) \, {\rm e}^{\frac {g_L \, t} C} = \frac C {g_L} \frac {g_L \, E_L +I} C \, {\rm e}^{\frac {g_L \, t} C} + \text{const}\\ ⟺ & \; V(t) \, {\rm e}^{\frac {g_L \, t} C} = \left(E_L + \frac I {g_L}\right) \, {\rm e}^{\frac {g_L \, t} C} + \text{const}\\ ⟺ & \; V(t) \, = \left(E_L + \frac I {g_L}\right) + \text{const} × {\rm e}^{\frac {g_L \, t} C}\\ \end{align*}\]And with $V(0) = E_L$:
\[∀t ≥0, \quad V(t) \, = E_L + \frac I {g_L}\left(1{\rm e}^{\frac {g_L \, t} C}\right)\]
The symmetry with respect to the $V = E_L$ between negative and positive current curves is clearly seen in this analytical expression.
As seen in Figure c. the numerical solution resulting from the Euler method and the exact one closely match.
So, on the whole:
 $V$ starts at $E_L$
 and then increases monotonously until $E_L + I/g_L$ (so the higher the input current, the higher the voltage, as examplified in Figure b)
 the characteristic time constant thereof is $C/g_L = 10 \text{ ms}$

the derivative of $V$:
\[\frac{ {\rm d}V}{ {\rm d}t} ≝ \frac{I}{C} \, {\rm e}^{\frac {g_L \, t} C}\]is positively proportional to $I$: the higher the input current, the steeper the increase of $V$.
Actionpotentialgenerating mechanism
To account for spike generation, let us assume that the neuron fires each time $V$ gets greater than a threshold $V_{th}$, after which the membrane voltage is reset to its initial value $E_L$. Let us set $V_{th} ≝ −63 \text{ mV}$:
For an input current of $1$ nA, we get $7$ spikes within the first $100$ ms. Let us investigate how the number of spikes varies with the input current:
For $I ≥ 0.75$ nA: The higher the input current, the steeper the increase of $V$ (as shown before) and hence the the faster $V$ reaches the threshold: so, as expected, we observe that the higher the input current, the higher the number of spikes within $100$ ms.
But for $I ≤ 0.7$ nA, no spike is fired. Indeed, the precise input current value at which the neuron start firing can be computed analytically:
\[\begin{align*} & \; ∃t; \; V(t) \, \overset{\text{eventually}}{>} V_{th}\\ ⟺ & \; \lim\limits_{t \to +∞} V(t) > V_{th} &&\text{ (as } V \text{ is nondecreasing)}\\ ⟺ & \; E_L + \frac I {g_L} > V_{th}\\ ⟺ & \; I > g_L (V_{th}  E_L) \end{align*}\]As it happens, with with $E_L ≝ 70 \text{ mV}$ and $V_{th} = 63 \text{ mV}$:
\[∃ t; \quad V(t) \, \overset{\text{eventually}}{>} V_{th} ⟺ I > 0.7 \text{ nA}\]
Let us plot a rastergram of the spike trains for various injected currents, so as to to see more easily how the number of spikes is impacted:
For $I ≥ 7.8$ nA, the number of spikes reaches its maximum value of $49$ spikes within $100$ ms.
To make it even more clear, let us plot the tuning curve of this neuron (i.e. the number of spikes within $100 \text{ ms}$ depending on the input current $I$):
As shown before, the current threshold at which the neuron starts firing is
\[I_{th} ≝ g_L (V_{th}  E_L)\]This makes sense intuitively:

the higher the difference $V_{th}  E_L$, the more the voltage $V$ has to increase to reach $V_{th}$ within $100 \text{ ms}$, and we saw previously that it requires a higher input current $I$

the higher the leak conductance $g_L$ of the membrane, the lower the voltage $V$ by Ohm’s law (as $V$ is inversely proportional to the conductance), and the higher the input current is required to be, to compensate
Refractory period
Whenever the neuron depolarize until the voltage hits the threshold: it fires, then the voltage is reinitialized at $E_L$, but it can be shown experimentally that the voltage doesn’t increase right away after having been reset: for a fixed time span (the refractory period), the voltage remains constant, then at the end of the refractory period it depolarizes again.
We can take into account this refractory period in our model:
On goes from $7$ spikes (before) to $5$ ones within the first $100$ ms, with a $5$ ms refractory period.
As expected: the higher the refractory period, the less spike the neuron fires within the same time span.
With a white noise term
A more realistic model would take into account part of the noise there is in reality. A more noisy voltage accross a neuron’s passive membrane could be given by:
\[C \frac{ {\rm d}V}{ {\rm d}t} = g_L (E_L  V(t))+ I + ση(t)\]where $σ$ is the noise magnitude and $η \sim 𝒩(0, 1)$. The choice of a Gaussian distribution is accounted for by the central limit theorem, since the noise can be thought of as a sum of many independent indentically distributed processes.
The Euler method yields:
\[C (V(t+ Δt)  V(t)) = g_L (E_L  V(t)) Δt + I Δt + σ \sqrt{Δt} \, η(t)\]Thus:
\[V(t+Δt) = \Big(1 \frac {Δt \cdot g_L} C\Big) V(t) + \frac {Δt \, (g_L \, E_L + I)} C + \frac {σ \sqrt{Δt} \, η(t)} C\]The $\sqrt{Δt}$ stems from the fact that
\[\sqrt{Δt} \, η(t) \sim 𝒩(0, Δt)\]so that the noise magnitude amounts to $Δt$ when adding a $Δt$ time step.
The white noise term makes it look more realistic:
2. IntegrateandFire neuron with timevarying input current to match experimental data
Finally, let us try to model the experminental spike trains generated by a given neuron in the primary somatosensory cortex of a monkey subject to a vibratory stimulus on the finger (as seen in the previous problem).
To do so, a constant input current will not be enough, as shown in Figure g.2.: the rastergram is not likely to match what is observed experimentally. And above all, the vibration frequency $f_1$ has not even been taken into account so far, which cannot but lead us to fail to model the experimental spike trains, as we showed that there is a linear relationship between the average spike counts/firing rates and the vibration frequencies.
Looking at the experimental spike trains, we observed in the previous problem that the spikes seemed to occur periodically within the simulation time span (between $200$ and $700$ ms): the higher the vibration frequency, the smaller the interspike period.
This suggests a sinusoidal input current whose frequency is a multiple of the vibration frequency, modulated by an envelope that decays out of the simulation period.
Gaussian envelope
With a Gaussian envelope, the timevarying input current
\[I(t) \propto \exp\left(\frac{(t  μ)^2}{2 s^2}\right) \cos(f_1 (tμ)/200)\]where $μ = 450 \text{ ms}, \; s = 150 \text{ ms}$, could look like this:
It leads to the following rastergram:
which is not satisfying, with regard to the experimental one.
Heaviside envelope
Another solution would be to go for a Heaviside envelope: the current is purely sinusoidal within the simulation period, it is equal to zero everywhere else.
The resulting timevarying input current
\[I(t) \propto \begin{cases} \cos(f_1 (t200)/200) \text{ if } t∈[200, 700] \text{ ms} \\ 0 \text{ else} \end{cases}\]could look like this:
The resulting rastergram:
resembles way more to the experimental one!
However, the average spike count/firing rate as a function of the vibration frequency is not satisfying: it is somewhat antilinear with respect to the vibration frequency!
This can be explained by the fact that even though the frequency of spikes increases with the vibration frequency $f_1$, the width of the “spikes fired in a row” is decreasing (since the sinusoidal current curve is becoming thinner and thinner, the current reaches high values during a shorter and shorter time as the frequency increases). To overcome that, we would have to ensure that the width of the current waves remains the same irrespective of the frequency: this can be achieved with square waves.
If the model somehow described what happens in reality, it would mean that the current injected though the membrane of the monkey’s primary somatosensory cortex neuron when a stimulus vibrates at a frequency $f_1$ on its fingertip is:
 null outside the simulation period
 sinusoidal within the simulation period
leading to a voltage accross the membrane reaching the threshold voltage periodically (the higher the vibration frequency, the smaller the period), which in turn induces an average firing rate depending linearly on $f_1$.
Leave a comment