# Colab setup ------------------
import os, sys, subprocess
if "google.colab" in sys.modules:
cmd = "pip install --upgrade iqplot watermark"
process = subprocess.Popen(cmd.split(), stdout=subprocess.PIPE, stderr=subprocess.PIPE)
stdout, stderr = process.communicate()
# ------------------------------import types
import numpy as np
import iqplot
import bokeh.io
import bokeh.plotting
bokeh.io.output_notebook()While not a story, per se, we can sometimes define probability distributions in terms of differential equations. Master equations and Fokker-Planck equations describe the dynamics, respectively, of probability mass functions and probability density functions over time. Langevin equations and stochastic differential equations describe quantities that may undergo random fluctuations.
As an example of sampling out of probability distributions described by differential equations, we will consider here a leaky integrate-and-fire (LIF) model which describes the membrane potential of a neuron. The membrane of a neuron is modeled as a capacitor which can accumulate charge and a resistor in parallel. For an input current \(I(t)\), then, the membrane potential \(V(t)\) obeys, via Kirchhoff’s law,
\[ \begin{aligned} I(t) = C\,\frac{\mathrm{d}V}{\mathrm{d}t} + \frac{V(t)}{R}, \end{aligned} \tag{15.1}\]
or, upon rearrangement,
\[ \begin{aligned} \frac{\mathrm{d}V}{\mathrm{d}t} = - \frac{V(t)}{RC} + \frac{I(t)}{C}. \end{aligned} \tag{15.2}\]
As the membrane potential grows with time, when it hits a threshhold \(V_\mathrm{thresh}\), an action potential is fired, after which the membrane potential rapidly relaxes back to a reset potential, \(V_\mathrm{reset}\). We define
\[ \begin{aligned} v(t) = \frac{V(t) - V_\mathrm{reset}}{V_\mathrm{thresh} - V_\mathrm{reset}} \end{aligned} \tag{15.3}\]
as a dimensionless membrane potential. Note that the reset potential corresponds to \(v = 0\) and the threshold corresponds to \(v = 1\). We further define a dimensionless current, \(i(t)\), such that \(I(t) = I_0\,i(t)\), where \(I_0\) is the amplitude of the input current. Using these expressions, we have
\[ \begin{aligned} \frac{\mathrm{d}v}{\mathrm{d}t} = - \frac{v(t)}{RC} - \frac{V_\mathrm{reset}}{(V_\mathrm{thresh} - V_\mathrm{reset})RC} + \frac{I_0\,i(t)}{(V_\mathrm{thresh} - V_\mathrm{reset})C}. \end{aligned} \tag{15.4}\]
There are two time scales we can identify. First, \(\tau_m = RC\) is the time scale of the potential drop of the membrane due to leakage. Second, \(\tau = (V_\mathrm{thresh} - V_\mathrm{reset})C / I_0\) is the amount of time it takes for the input current to bring the membrane up to the threshold potential from the reset potential. With these definitions, and also defining \(\alpha = V_\mathrm{reset}/(V_\mathrm{thresh} - V_\mathrm{reset})\), we have
\[ \begin{aligned} \frac{\mathrm{d}v}{\mathrm{d}t} = - \frac{v(t) + \alpha}{\tau_m} + \frac{i(t)}{\tau}. \end{aligned} \tag{15.5}\]
This is a deterministic differential equation. We may have noise in the input current, which we include in the equation as \(\sigma\,\xi(t)\), giving us a Langevin equation
\[ \begin{aligned} \frac{\mathrm{d}v}{\mathrm{d}t} = - \frac{v(t) + \alpha}{\tau_m} + \frac{i(t)}{\tau} + \sigma\,\xi(t). \end{aligned} \tag{15.6}\]
The function \(i(t)\) can be arbitrarily complicated. In many models, it is itself stochastic, modeling input currents from neighboring neurons. Similarly, the noise function \(\xi(t)\) may also take a variety of forms. We will focus on Gaussian white noise, which means that the stochastic portion of the dynamics of the voltage is fluctuating according to a Normal distribution with zero mean. (More precisely, \(\xi(t)\) is defined by a Wiener process, which means, among other things, that \(\langle \xi(t)\rangle = 0\) and \(\langle \xi(t)\,\xi(t')\rangle = \delta(t-t')\).)
We will solve the above stochastic LIF Langevin equation using Euler-Maruyama integration. The idea is that we integrate exactly like Euler’s method with time step size \(dt\), but we add a random number drawn from \(\text{Norm}(0, dt)\) (which is then multiplied by \(\sigma\), the scale of the noise) at every step. Note that this is equivalent to drawing out of a standard normal (with mean zero and variance of one) and multiplying the result by \(\sigma\,dt\). In this case, every time the dimensionless voltage reaches the spike threshold (\(v = 1\)), we record a spike and then immediately send the voltage back down to the reset voltage (\(v = 0\)). This is implemented in the function below.
# Instantiate random number generator
rng = np.random.default_rng()
def integrate_noisy_lif(
v0=0.0,
dt=0.001,
T=1000.0,
tau_m=1000.0,
tau=30.0,
alpha=1.0,
sigma=1.0,
ifun=1.0,
ifun_args=(),
return_traj=False,
):
"""Integrate the stochastic LIF model and return spike times.
Parameters
----------
v0 : float, default 0.0
Dimensionless voltage at time t = 0.
dt : float, default 0.001
Time step interval.
T : float, default 1000.0
Integration proceeds from t = 0 to t = T.
tau_m : float, default 1000.0
Time scale of membrane potential relaxation due to leakage.
tau : float, default 30.0
Time scale for membrane potential to rise to threshold.
alpha : float, default 1.0
alpha parameter from LIF model.
sigma : float, default 1.0
Magnitude of noise.
ifun : function or float, default 1.0
Function with call signature `ifun(t, *ifun_args)` that gives
the current as a function of time. If specified as a float,
it is a constant function.
ifun_args : tuple, default ()
Arguments to pass into `ifun`. Ignored if `ifun` is a float.
return_traj : bool, default False
If True, also return (t, v), the time points and corresponding
dimensionless voltage, as a tuple in addition to spike times.
Returns
-------
spikes : Numpy array
Times of spikes.
traj_tuple : tuple of Numpy arrays (optional)
If `return_traj` is True, a tuple (t, v) of the time and
corresponding dimensionless voltage.
Notes
-----
.. The parameters `dt`, `T`, `tau_m`, and `tau` all must be in the
same units. These are the same units of the time array returned
if `traj_tuple` is True.
.. The default parameter values of `tau` and `tau_m` are such that
the time units are in milliseconds and are in the limit of very
slow leakage and a "typical" firing rate of about 30 spikes/s.
"""
# Set up ifun, it not given as a function
if type(ifun) != types.FunctionType:
ret_val = float(ifun)
ifun = lambda t: ret_val
# Set up time points
t = np.arange(0, T, dt)
# Output array
v = np.empty_like(t)
v[0] = v0
# Initialize list to store spikes
spikes = []
for i in range(1, len(v)):
# Euler-Maruyama step
v[i] = v[i - 1] + dt * (
-(v[i - 1] + alpha) / tau_m
+ ifun(t[i - 1], *ifun_args) / tau
+ sigma * rng.standard_normal()
)
# Record spike and reset voltage if we passed threshold
if v[i] >= 1.0:
spikes.append(t[i])
v[i] = 0.0
if return_traj:
return np.array(spikes), (t, v)
else:
return np.array(spikes)Let’s give this function a while to sample some spike times!
spikes, (t, v) = integrate_noisy_lif(return_traj=True)We can plot the trajectory.
p = bokeh.plotting.figure(
frame_width=500,
frame_height=200,
x_axis_label='t (ms)',
y_axis_label='v(t)'
)
# Thin the trajectory a bit when plotting
p.line(t[::10], v[::10])
bokeh.io.show(p)We can see what happens if we increase the leakage, making the two time scales \(\tau_m\) and \(\tau\) comparable. We will try \(\tau_m = 70\) ms and \(\tau = 30\) ms.
spikes, (t, v) = integrate_noisy_lif(tau_m=70.0, return_traj=True)
p = bokeh.plotting.figure(
frame_width=500,
frame_height=200,
x_axis_label='t (ms)',
y_axis_label='v(t)'
)
# Thin the trajectory a bit when plotting
p.line(t[::10], v[::10])
bokeh.io.show(p)We get far less frequent spikes. This is because the leakage is causing the potential to drop as it gets large.
Ultimately our goal is to sample spike times (actually, interspike intervals, or ISIs), so we can run a longer trajectory to get plenty of samples.
spikes = integrate_noisy_lif(T=10_000)
# Look as a strip plot
bokeh.io.show(
iqplot.strip(
spikes,
marker="dash",
marker_kwargs=dict(line_width=0.5),
x_axis_label="time (ms)",
frame_height=100,
frame_width=700,
)
)We can compute interspike intervals by computing the difference between successive spike times and can plot the ECDF.
bokeh.io.show(
iqplot.ecdf(np.diff(spikes), x_axis_label='ISI (ms)')
)As we will see in the exercises, with the above choice of parameters, the ISI for the noisy leaky integrate-and-fire model is approximately Inverse Gaussian distributed. However, ISIs are often Exponentially distributed, which means that the spiking is a Poisson process. As described by Stevens and Zador, the noise leaky integrate-and-fire model can give approximately Exponentially distributed interspike intervals when the input current itself is also given by Gaussian white noise with a positive mean. They further showed that this works when spiking is slow, such that the characteristic interspike interval is much greater than \(\tau_m\).
As an example, consider the case where \(\tau_m = 70\) ms and \(\tau = 40\) ms and the current function has a mean of one (in dimensionless units) and a standard deviation of 0.3. We can code this up below.
def ifun(t, sigma_I=0.3):
return rng.normal(1.0, sigma_I)
spikes = integrate_noisy_lif(tau_m=70, tau=40, ifun=ifun, T=100_000)To demonstrate that the ISIs are approximately Exponentially distributed, we plot the empirical complementary cumulative distribution function (ECCDF), which is \(1 - \hat{F}(t_\mathrm{ISI})\), where \(\hat{F}(t_\mathrm{ISI})\) is the ECDF. Importantly, a linear tail on a semilog plot of the ECCDF indicates an Exponential distribution, since the CDF for the Exponential distribution is \(F(t_\mathrm{ISI}) = 1 - \mathrm{e}^{-\beta t_\mathrm{ISI}}\) such that the ECCDF is \(\mathrm{e}^{-\beta t_\mathrm{ISI}}\).
# Plot ECCDF
p = iqplot.ecdf(np.diff(spikes), x_axis_label='ISI (ms)', complementary=True, y_axis_type='log')
# Annotate with a slope 1 line
t_shift = 110
beta = 1 / np.mean(np.diff(spikes) - t_shift)
x = np.array([110, 1500])
y = np.exp(-beta * (x - t_shift))
p.line(x, y, color='tomato', line_width=2)
bokeh.io.show(p)So, for ISIs large compared to \(\tau_m\), the distribution of ISIs is approximately Exponential. In practice, this is like an Exponential distribution with a time lag. That is, spikes arrive a Poisson processes, but only after a refractory period after the previous spike. We will investigate a phenomenological model for Poissonian-with-refractory-period spikes in the exercises.
Spiking may also appear Poissonian if the stimuli (the input current) arrive in a Poissonian manner and the LIF circuit responds quickly (\(\tau\) is small). This is kind of obvious; a circuit that responds close to instantaneously to a Poissionian stimulus will itself behave in a Poissonian manner.
To demonstrate this, we can put in current as Poissonian-distributed spikes with a Gaussian shape.
# Arrival of stimulating current as a Poisson process
stimulus_times = np.cumsum(rng.exponential(30, size=1000))
def ifun(t, stim_times, I_peak=1.0, stim_width=1.0):
"""Stimuli come as Gaussian-like peaks"""
return np.sum(I_peak * np.exp(-(t-stimulus_times)**2 / 2 / stim_width**2))
spikes = integrate_noisy_lif(tau=3, ifun=ifun, ifun_args=(stimulus_times,), T=30_000)We can compare the ISIs to what we would expect from a Poisson process (exponentially distributed ISI).
p = iqplot.ecdf(np.diff(spikes), x_axis_label='ISI (ms)')
x = np.linspace(0, 275, 300)
y = 1 - np.exp(-x / np.mean(np.diff(spikes)))
p.line(x, y, line_width=2, color='tomato')
bokeh.io.show(p)Indeed, we have Exponentially distributed ISIs.
%load_ext watermark
%watermark -v -p numpy,iqplot,bokeh,jupyterlabPython implementation: CPython
Python version : 3.12.9
IPython version : 9.1.0
numpy : 2.1.3
iqplot : 0.3.7
bokeh : 3.6.2
jupyterlab: 4.4.2