Example: Tsodyks_Uziel_Markram_2000
Note
You can launch an interactive, editable version of this example without installing any local files using the Binder service (although note that at some times this may be slow or fail to open):
Fig. 1 from:
Synchrony Generation in Recurrent Networks with Frequency-Dependent Synapses The Journal of Neuroscience, 2000, Vol. 20 RC50
Implementation partially based on nest-2.0.0/examples/nest/tsodyks_shortterm_bursts.sli by Moritz Helias, 2006.
Sebastian Schmitt, 2022
import numpy as np
# set seed for reproducible figures
np.random.seed(5)
# for truncated normal
import scipy
from scipy import stats
import matplotlib.pyplot as plt
from brian2 import (
NeuronGroup,
Synapses,
SpikeGeneratorGroup,
SpikeMonitor,
StateMonitor,
)
from brian2 import ms, mV
from brian2 import run, defaultclock
def truncated_normal(loc, scale, bounds, size):
"""Normal distribution truncated within bounds
loc -- mean (“centre”) of the distribution
scale -- standard deviation (spread or “width”) of the distribution
bounds -- list of min and maximum
size -- number of samples
"""
bounds = np.array([bounds] * size)
s = scipy.stats.truncnorm.rvs(
(bounds[:, 0] - loc) / scale, (bounds[:, 1] - loc) / scale, loc=loc, scale=scale
)
return s
def get_population(name, N, tau_refrac):
"""Get population of neurons
name -- name of population
N -- number of neurons
tau_refrac -- refractory period
"""
neurons = NeuronGroup(
N,
"""
tau_mem : second
tau_refrac : second
v_reset : volt
v_thresh : volt
I_syn_ee_synapses : volt
I_syn_ei_synapses : volt
I_syn_ie_synapses : volt
I_syn_ii_synapses : volt
I_b : volt
dv/dt = -v/tau_mem + (I_syn_ee_synapses +
I_syn_ei_synapses +
I_syn_ie_synapses +
I_syn_ii_synapses)/tau_mem
+ I_b/tau_mem : volt (unless refractory)
""",
threshold="v>v_thresh",
reset="v=v_reset",
refractory=tau_refrac,
method="exact",
name=name,
)
v_thresh = 15 * mV
v_reset = 13.5 * mV
neurons.tau_mem = 30 * ms
neurons.v_thresh = v_thresh
neurons.v_reset = v_reset
# paper gives range of 0.05 mV but population bursts are not visible with that value
# -> increased to 1 mV range
neurons.I_b = (
np.random.uniform(v_thresh / mV - 0.5, v_thresh / mV + 0.5, size=N) * mV
)
return neurons
def get_synapses(name, source, target, tau_I, A, U, tau_rec, tau_facil=None):
"""Construct connections and retrieve synapses
name -- name of synapses
source -- source of connections
target -- target of connections
tau_I -- inactivation time constant
A -- absolute synaptic strength
U -- utilization of synaptic efficacy
tau_rec -- recovery time constant
tau_facil -- facilitation time constant (optional)
"""
synapses_eqs = """
A : volt
U : 1
tau_I : second
tau_rec : second
dx/dt = z/tau_rec : 1 (clock-driven) # recovered
dy/dt = -y/tau_I : 1 (clock-driven) # active
z = 1 - x - y : 1 # inactive
I_syn_{}_post = A*y : volt (summed)
""".format(
name
)
if tau_facil:
synapses_eqs += """
du/dt = -u/tau_facil : 1 (clock-driven)
tau_facil : second
"""
synapses_action = """
u += U*(1-u)
y += u*x # important: update y first
x += -u*x
"""
else:
synapses_action = """
y += U*x # important: update y first
x += -U*x
"""
synapses = Synapses(
source,
target,
model=synapses_eqs,
on_pre=synapses_action,
method="exact",
name=name,
)
synapses.connect(p=0.1)
N_syn = len(synapses)
synapses.tau_I = tau_I
A_min = min(0.2 * A, 2 * A)
A_max = max(0.2 * A, 2 * A)
synapses.A = (
truncated_normal(
A / mV, 0.5 * abs(A / mV), [A_min / mV, A_max / mV], size=N_syn
) * mV
)
assert not any(synapses.A < A_min)
assert not any(synapses.A > A_max)
U_mean, U_min, U_max = U
synapses.U = truncated_normal(U_mean, 0.5 * U_mean, [U_min, U_max], size=N_syn)
assert not any(synapses.U <= U_min)
assert not any(synapses.U > U_max)
tau_min = 5
synapses.tau_rec = (
truncated_normal(
tau_rec / ms, 0.5 * tau_rec / ms, [tau_min, np.inf], size=N_syn
) * ms
)
assert not any(synapses.tau_rec / ms <= tau_min)
if tau_facil:
synapses.tau_facil = (
truncated_normal(
tau_facil / ms, 0.5 * tau_facil / ms, [tau_min, np.inf], size=N_syn
) * ms
)
assert not any(synapses.tau_facil / ms <= tau_min)
# start fully recovered
synapses.x = 1
return synapses
# configure neuron populations
exc_neurons = get_population("exc_neurons", N=400, tau_refrac=3 * ms)
inh_neurons = get_population("inh_neurons", N=100, tau_refrac=2 * ms)
# configure synapses
ee_synapses = get_synapses(
"ee_synapses",
exc_neurons,
exc_neurons,
tau_I=3 * ms,
A=1.8 * mV,
U=[0.5, 0.1, 0.9],
tau_rec=800 * ms,
)
ei_synapses = get_synapses(
"ei_synapses",
exc_neurons,
inh_neurons,
tau_I=3 * ms,
A=7.2 * mV,
U=[0.04, 0.001, 0.07],
tau_rec=100 * ms,
tau_facil=1000 * ms,
)
ie_synapses = get_synapses(
"ie_synapses",
inh_neurons,
exc_neurons,
tau_I=3 * ms,
A=-5.4 * mV,
U=[0.5, 0.1, 0.9],
tau_rec=800 * ms,
)
ii_synapses = get_synapses(
"ii_synapses",
inh_neurons,
inh_neurons,
tau_I=3 * ms,
A=-7.2 * mV,
U=[0.04, 0.001, 0.07],
tau_rec=100 * ms,
tau_facil=1000 * ms,
)
# run for burnin time to settle network activity
defaultclock.dt = 1 * ms
burnin = 900
run(burnin * ms)
# record from now on
spike_monitor_exc = SpikeMonitor(exc_neurons)
spike_monitor_inh = SpikeMonitor(inh_neurons)
state_monitor_ee = StateMonitor(ee_synapses, ["x"], record=True)
duration = 4200
run(duration * ms, report="text")
# plots
fig, axes = plt.subplots(3, figsize=(6, 8), sharex=True)
# raster plot
axes[0].plot(spike_monitor_exc.t / ms, spike_monitor_exc.i, ".k", ms=1)
axes[0].plot(spike_monitor_inh.t / ms, spike_monitor_inh.i + len(exc_neurons), ".k", ms=1)
axes[0].set_ylabel("Neuron No.")
axes[0].set_ylim(0, len(exc_neurons) + len(inh_neurons))
# network activity
net_activity = np.histogram(
np.concatenate(
list(spike_monitor_exc.spike_trains().values())
+ list(spike_monitor_inh.spike_trains().values())
) / ms,
bins=np.arange(burnin, duration + burnin, 1))[0] / (len(exc_neurons) + len(inh_neurons))
axes[1].plot(np.arange(0, len(net_activity)) + burnin, net_activity, "k")
net_activity_min = 0
net_activity_max = 0.2
axes[1].set_ylim(net_activity_min, net_activity_max)
axes[1].set_ylabel("Net activity")
# network activity inset
axins = axes[1].inset_axes([0.05, 0.35, 0.2, 0.6])
axins.plot(np.arange(0, len(net_activity)) + burnin, net_activity, "k")
inset_min = 1220
inset_max = 1260
axins.set_xlim(inset_min + burnin, inset_max + burnin)
axins.set_ylim(net_activity_min, net_activity_max)
axins.set_xticks([inset_min + burnin, inset_max + burnin])
axins.set_xticklabels([inset_min, inset_max])
axins.set_yticks([])
# recovered synaptic partition
axes[2].plot(
state_monitor_ee.t / ms, np.mean(state_monitor_ee.x, axis=0), "k", label="x"
)
axes[2].set_ylim(0.2, 0.6)
axes[2].set_xlabel("Time (msec)")
axes[2].set_ylabel("Recov excit")
axes[2].set_xlim(burnin, duration + burnin)
xtickstep = 1000
axes[2].set_xticks(np.arange(burnin, duration + burnin, xtickstep))
axes[2].set_xticklabels(map(str, range(0, duration, xtickstep)))
axes[0].xaxis.set_tick_params(which="both", labelbottom=True)
axes[1].xaxis.set_tick_params(which="both", labelbottom=True)
plt.show()
