Note. The following notebook contains code in addition to text and figures. By default, the code has been hidden. You can click the icon that looks like an eye in the toolbar above to show the code. To run the code, click the cell menu, then "run all".
In [ ]:
# Import packages, set preferences, etc.
%matplotlib inline
from brian2 import *
import ipywidgets as ipw
from numpy.random import poisson
import warnings
warnings.filterwarnings("ignore")
prefs.codegen.target = 'cython'
defaultclock.dt = 0.05*ms
In [ ]:
%%html
<!-- hack to improve styling of ipywidgets sliders -->
<style type="text/css">
.widget-label {
min-width: 35ex;
max-width: 35ex;
}
.widget-hslider {
width: 100%;
}
</style>
We start from the following model of a chopper cell in the ventral cochlear nucleus. The cell is a leaky integrate and fire neuron, defined by the differential equation $$\tau\frac{\mathrm{d}v}{\mathrm{d}t}=-v.$$ If $v>1$ it fires a spike and instantly resets to $v\leftarrow 0$ and is held there for time $t_\mathrm{ref}$. It receives $N$ excitatory and $N$ inhibitory synapses from neurons that fire Poisson spikes at rates $\rho_E$ and $\rho_I$ respectively. The excitatory synapses have weight $w$ and the inhibitory $-w$, so an incoming spike causes an instantaneous change $v\leftarrow v\pm w$. We will sometimes write the inhibitory fraction $\alpha=\rho_I/\rho_E$ and reparametrise the pair $(\rho_E, \rho_I)$ as $(\rho, \alpha)$ where $\rho_E=\rho$ and $\rho_I=\alpha\rho$.
The following interactive figure allows you to visualise the behaviour of this model for different values of the parameters $\tau$, $N$, $w$, $\rho$, $\alpha$, $t_\mathrm{ref}$.
In [ ]:
def basic_model(N=50, tau_ms=10, w=0.02, rho_Hz=200, alpha=0.0, t_ref_ms=0.1):
# Set parameters
repeats = 100
duration = 100*ms
tau = tau_ms*ms
rho = rho_Hz*Hz
t_ref = t_ref_ms*ms
# Define and run the model
eqs = '''
dv/dt = -v/tau : 1 (unless refractory)
'''
G = NeuronGroup(repeats, eqs, threshold='v>1', reset='v=0',
refractory=t_ref, method='euler')
Pe = PoissonInput(G, 'v', N, rho, weight=w)
Pi = PoissonInput(G, 'v', N, rho*alpha, weight=-w)
spikemon = SpikeMonitor(G)
statemon = StateMonitor(G, 'v', record=[0])
run(duration)
# Compute ISI histograms
isi = []
for train in spikemon.spike_trains().values():
train.sort()
isi.append(diff(train))
isi = hstack(isi)
cv = std(isi)/mean(isi)
# Plot results
figure(figsize=(10, 2.5))
subplot(131)
plot(spikemon.t/ms, spikemon.i, '.k')
xlabel('Time (ms)')
ylabel('Repeat number')
title('Spike raster plot')
xlim(0, duration/ms)
ylim(0, repeats)
subplot(132)
plot(statemon.t/ms, statemon.v.T, '-k')
xlabel('Time (ms)')
ylabel('v')
title('Membrane potential trace')
xlim(0, duration/ms)
ylim(-0.2, 1.2)
axhline(0, ls=':', c='r')
axhline(1, ls=':', c='g')
subplot(133)
hist(isi/ms, fc='k', bins=arange(60)*0.5)
yticks([])
ylabel('Frequency')
xlabel('ISI (ms)')
title('Interspike interval histogram')
#title('CV = %.2f' % cv)
text(0.95, 0.9, 'CV = %.2f' % cv, ha='right', va='top',
bbox=dict(facecolor='white'),
transform=gca().transAxes)
tight_layout()
display(ipw.interact(basic_model,
N=ipw.IntSlider(
min=1, max=100, step=1, value=50,
continuous_update=False,
description="Number of inputs $N$"),
tau_ms=ipw.FloatSlider(
min=0.1, max=20.0, step=0.1, value=10.0,
continuous_update=False,
description=r"Membrane time constant $\tau$ (ms)"),
rho_Hz=ipw.FloatSlider(
min=10, max=500, step=10, value=200.0,
continuous_update=False,
description=r"Excitatory input firing rate $\rho$ (sp/s)"),
w=ipw.FloatSlider(
min=0, max=1, step=0.01, value=0.02,
continuous_update=False,
description=r"Synaptic weight $w$"),
alpha=ipw.FloatSlider(
min=0, max=1.5, step=0.05, value=0,
continuous_update=False,
description=r"Inhibitory fraction $\alpha$"),
t_ref_ms=ipw.FloatSlider(
min=0, max=5, step=0.05, value=0.1,
continuous_update=False,
description=r"Refractory period $t_\mathrm{ref}$ (ms)"),
));
We now simplify the model using the diffusion approximation of Fourcaud and Brunel (2002). The equations above can be approximated by a leaky integrate and fire neuron with the stochastic differential equation $$\tau\frac{\mathrm{d}v}{\mathrm{d}t}=-v+\mu+\sigma\sqrt{\tau}\xi,$$ where $\xi$ is a stochastic differential and the constants are computed as $$\begin{align} \mu &= wN\tau(\rho_E-\rho_I)=wN\tau\rho(1-\alpha)\\ \sigma^2 &= w^2N\tau(\rho_E+\rho_I)=w^2N\tau\rho(1+\alpha). \end{align}$$
This approximation eliminates the need for the input neurons and synaptic connections. You can see in the interactive figure below how well this approximation holds. As long as the number of inputs $N$ and the firing rate $\rho$ are reasonably large, the two curves are statistically very similar (although of course since each trajectory is random, they are not close point-by-point and cannot be). Note that in the figure below, the spiking behaviour when $v>1$ has been removed to see the qualitative behaviour of the stochastic differential equation more clearly.
In [ ]:
def diffusion_panel(N=50, tau_ms=10, w=0.02, rho_Hz=200, alpha=0.0):
# Parameters
tau = tau_ms*ms
rho = rho_Hz*Hz
duration = 200*ms
# diffusion approximation parameters
mu = w*N*tau*rho*(1-alpha)
sigma = w*sqrt(N*tau*rho*(1+alpha))
# Define and run the model
eqs = '''
dv/dt = -v/tau : 1
du/dt = (mu-u)/tau+sigma*xi*tau**-0.5 : 1
'''
G = NeuronGroup(1, eqs, method='euler')
Pe = PoissonInput(G, 'v', N, rho, weight=w)
Pi = PoissonInput(G, 'v', N, rho*alpha, weight=-w)
M = StateMonitor(G, ['v', 'u'], record=[0])
run(duration)
# Plot
plot(M.t/ms, M.v[0], label="Original")
plot(M.t/ms, M.u[0], label="Diffusion approximation")
yticks([])
xlabel('Time (ms)')
ylabel('Membrane potential')
legend(loc='lower right')
display(ipw.interact(diffusion_panel,
N=ipw.IntSlider(
min=1, max=100, step=1, value=50,
continuous_update=False,
description=r"Number of inputs $N$:"),
tau_ms=ipw.FloatSlider(
min=0.1, max=20.0, step=0.1, value=10.0,
continuous_update=False,
description=r"Membrane time constant $\tau$ (ms):"),
rho_Hz=ipw.FloatSlider(
min=10, max=500, step=10, value=200.0,
continuous_update=False,
description=r"Excitatory input firing rate $\rho$ (sp/s):"),
w=ipw.FloatSlider(
min=0, max=1, step=0.01, value=0.02,
continuous_update=False,
description=r"Synaptic weight $w$:"),
alpha=ipw.FloatSlider(
min=0, max=1.5, step=0.05, value=0,
continuous_update=False,
description=r"Inhibitory fraction $\alpha$:"),
));
In the diffusion approximation, all the parameters $N$, $\rho_E$, $\rho_I$, $w$ reduce to just $\mu$ and $\sigma$ meaning that they are not independent. It is possible to directly set these values and this gives us a model with four independent parameters $\mu$, $\sigma$, $\tau$, and $t_\mathrm{ref}$. We call this the reduced model and you can explore its behaviour with the interactive figure below. Note that the larger you make $\sigma$ compared to $\mu$, the higher the CV (typically). We will explore this more in the next notebook, Maps.
In [ ]:
def reduced_model(mu=1.5, sigma=0.5, tau_ms=10, t_ref_ms=0.1):
# Set parameters
repeats = 100
duration = 100*ms
tau = tau_ms*ms
t_ref = t_ref_ms*ms
# Define and run the model
eqs = '''
dv/dt = (mu-v)/tau+sigma*xi*tau**-0.5 : 1
'''
G = NeuronGroup(repeats, eqs, threshold='v>1', reset='v=0',
refractory=t_ref, method='euler')
spikemon = SpikeMonitor(G)
statemon = StateMonitor(G, 'v', record=[0])
run(duration)
# Compute ISI histograms
isi = []
for train in spikemon.spike_trains().values():
train.sort()
isi.append(diff(train))
isi = hstack(isi)
cv = std(isi)/mean(isi)
# Plot results
figure(figsize=(10, 2.5))
subplot(131)
plot(spikemon.t/ms, spikemon.i, '.k')
xlabel('Time (ms)')
ylabel('Repeat number')
title('Spike raster plot')
xlim(0, duration/ms)
ylim(0, repeats)
subplot(132)
plot(statemon.t/ms, statemon.v.T, '-k')
xlabel('Time (ms)')
ylabel('v')
title('Membrane potential trace')
xlim(0, duration/ms)
ylim(-0.2, 1.2)
axhline(0, ls=':', c='r')
axhline(1, ls=':', c='g')
subplot(133)
hist(isi/ms, fc='k', bins=arange(60)*0.5)
yticks([])
ylabel('Frequency')
xlabel('ISI (ms)')
title('Interspike interval histogram')
#title('CV = %.2f' % cv)
text(0.95, 0.9, 'CV = %.2f' % cv, ha='right', va='top',
bbox=dict(facecolor='white'),
transform=gca().transAxes)
tight_layout()
display(ipw.interact(reduced_model,
tau_ms=ipw.FloatSlider(
min=0.1, max=20.0, step=0.1, value=10.0,
continuous_update=False,
description=r"Membrane time constant $\tau$ (ms):"),
t_ref_ms=ipw.FloatSlider(
min=0, max=5, step=0.05, value=0.1,
continuous_update=False,
description=r"Refractory period $t_\mathrm{ref}$ (ms):"),
mu=ipw.FloatSlider(
min=0, max=5, step=0.05, value=1.5,
continuous_update=False,
description=r"Mean current $\mu$:"),
sigma=ipw.FloatSlider(
min=0, max=5, step=0.05, value=0.5,
continuous_update=False,
description=r"Standard deviation of current $\sigma$:"),
));