1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
|
#!/usr/bin/env python3
"""
Fig. 8 from:
Brunel, N. Dynamics of Sparsely Connected Networks of Excitatory and
Inhibitory Spiking Neurons. J Comput Neurosci 8, 183–208
(2000). https://doi.org/10.1023/A:1008925309027
Inspired by http://neuronaldynamics.epfl.ch
Sebastian Schmitt, 2022
"""
import random
from brian2 import *
import matplotlib.pyplot as plt
def sim(g, nu_ext_over_nu_thr, sim_time, ax_spikes, ax_rates, rate_tick_step):
"""
g -- relative inhibitory to excitatory synaptic strength
nu_ext_over_nu_thr -- ratio of external stimulus rate to threshold rate
sim_time -- simulation time
ax_spikes -- matplotlib axes to plot spikes on
ax_rates -- matplotlib axes to plot rates on
rate_tick_step -- step size for rate axis ticks
"""
# network parameters
N_E = 10000
gamma = 0.25
N_I = round(gamma * N_E)
N = N_E + N_I
epsilon = 0.1
C_E = epsilon * N_E
C_ext = C_E
# neuron parameters
tau = 20 * ms
theta = 20 * mV
V_r = 10 * mV
tau_rp = 2 * ms
# synapse parameters
J = 0.1 * mV
D = 1.5 * ms
# external stimulus
nu_thr = theta / (J * C_E * tau)
defaultclock.dt = 0.1 * ms
neurons = NeuronGroup(N,
"""
dv/dt = -v/tau : volt (unless refractory)
""",
threshold="v > theta",
reset="v = V_r",
refractory=tau_rp,
method="exact",
)
excitatory_neurons = neurons[:N_E]
inhibitory_neurons = neurons[N_E:]
exc_synapses = Synapses(excitatory_neurons, target=neurons, on_pre="v += J", delay=D)
exc_synapses.connect(p=epsilon)
inhib_synapses = Synapses(inhibitory_neurons, target=neurons, on_pre="v += -g*J", delay=D)
inhib_synapses.connect(p=epsilon)
nu_ext = nu_ext_over_nu_thr * nu_thr
external_poisson_input = PoissonInput(
target=neurons, target_var="v", N=C_ext, rate=nu_ext, weight=J
)
rate_monitor = PopulationRateMonitor(neurons)
# record from the first 50 excitatory neurons
spike_monitor = SpikeMonitor(neurons[:50])
run(sim_time, report='text')
ax_spikes.plot(spike_monitor.t / ms, spike_monitor.i, "|")
ax_rates.plot(rate_monitor.t / ms, rate_monitor.rate / Hz)
ax_spikes.set_yticks([])
ax_spikes.set_xlim(*params["t_range"])
ax_rates.set_xlim(*params["t_range"])
ax_rates.set_ylim(*params["rate_range"])
ax_rates.set_xlabel("t [ms]")
ax_rates.set_yticks(
np.arange(
params["rate_range"][0], params["rate_range"][1] + rate_tick_step, rate_tick_step
)
)
plt.subplots_adjust(hspace=0)
parameters = {
"A": {
"g": 3,
"nu_ext_over_nu_thr": 2,
"t_range": [500, 600],
"rate_range": [0, 6000],
"rate_tick_step": 1000,
},
"B": {
"g": 6,
"nu_ext_over_nu_thr": 4,
"t_range": [1000, 1200],
"rate_range": [0, 400],
"rate_tick_step": 100,
},
"C": {
"g": 5,
"nu_ext_over_nu_thr": 2,
"t_range": [1000, 1200],
"rate_range": [0, 200],
"rate_tick_step": 50,
},
"D": {
"g": 4.5,
"nu_ext_over_nu_thr": 0.9,
"t_range": [1000, 1200],
"rate_range": [0, 250],
"rate_tick_step": 50,
},
}
for panel, params in parameters.items():
fig = plt.figure(figsize=(4, 5))
fig.suptitle(panel)
gs = fig.add_gridspec(ncols=1, nrows=2, height_ratios=[4, 1])
ax_spikes, ax_rates = gs.subplots(sharex="col")
sim(
params["g"],
params["nu_ext_over_nu_thr"],
params["t_range"][1] * ms,
ax_spikes,
ax_rates,
params["rate_tick_step"],
)
plt.show()
|
