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
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
|
# coding=utf-8
"""
Modeling neuron-glia interactions with the Brian 2 simulator
Marcel Stimberg, Dan F. M. Goodman, Romain Brette, Maurizio De PittÃ
bioRxiv 198366; doi: https://doi.org/10.1101/198366
Figure 1: Modeling of neurons and synapses.
Randomly connected networks with conductance-based synapses (COBA; see Brunel,
2000). Synapses exhibit short-time plasticity (Tsodyks, 2005; Tsodyks et al.,
1998).
"""
from brian2 import *
import sympy
import plot_utils as pu
seed(11922) # to get identical figures for repeated runs
################################################################################
# Model parameters
################################################################################
### General parameters
duration = 1.0*second # Total simulation time
sim_dt = 0.1*ms # Integrator/sampling step
N_e = 3200 # Number of excitatory neurons
N_i = 800 # Number of inhibitory neurons
### Neuron parameters
E_l = -60*mV # Leak reversal potential
g_l = 9.99*nS # Leak conductance
E_e = 0*mV # Excitatory synaptic reversal potential
E_i = -80*mV # Inhibitory synaptic reversal potential
C_m = 198*pF # Membrane capacitance
tau_e = 5*ms # Excitatory synaptic time constant
tau_i = 10*ms # Inhibitory synaptic time constant
tau_r = 5*ms # Refractory period
I_ex = 150*pA # External current
V_th = -50*mV # Firing threshold
V_r = E_l # Reset potential
### Synapse parameters
w_e = 0.05*nS # Excitatory synaptic conductance
w_i = 1.0*nS # Inhibitory synaptic conductance
U_0 = 0.6 # Synaptic release probability at rest
Omega_d = 2.0/second # Synaptic depression rate
Omega_f = 3.33/second # Synaptic facilitation rate
################################################################################
# Model definition
################################################################################
# Set the integration time (in this case not strictly necessary, since we are
# using the default value)
defaultclock.dt = sim_dt
### Neurons
neuron_eqs = '''
dv/dt = (g_l*(E_l-v) + g_e*(E_e-v) + g_i*(E_i-v) +
I_ex)/C_m : volt (unless refractory)
dg_e/dt = -g_e/tau_e : siemens # post-synaptic exc. conductance
dg_i/dt = -g_i/tau_i : siemens # post-synaptic inh. conductance
'''
neurons = NeuronGroup(N_e + N_i, model=neuron_eqs,
threshold='v>V_th', reset='v=V_r',
refractory='tau_r', method='euler')
# Random initial membrane potential values and conductances
neurons.v = 'E_l + rand()*(V_th-E_l)'
neurons.g_e = 'rand()*w_e'
neurons.g_i = 'rand()*w_i'
exc_neurons = neurons[:N_e]
inh_neurons = neurons[N_e:]
### Synapses
synapses_eqs = '''
# Usage of releasable neurotransmitter per single action potential:
du_S/dt = -Omega_f * u_S : 1 (event-driven)
# Fraction of synaptic neurotransmitter resources available:
dx_S/dt = Omega_d *(1 - x_S) : 1 (event-driven)
'''
synapses_action = '''
u_S += U_0 * (1 - u_S)
r_S = u_S * x_S
x_S -= r_S
'''
exc_syn = Synapses(exc_neurons, neurons, model=synapses_eqs,
on_pre=synapses_action+'g_e_post += w_e*r_S')
inh_syn = Synapses(inh_neurons, neurons, model=synapses_eqs,
on_pre=synapses_action+'g_i_post += w_i*r_S')
exc_syn.connect(p=0.05)
inh_syn.connect(p=0.2)
# Start from "resting" condition: all synapses have fully-replenished
# neurotransmitter resources
exc_syn.x_S = 1
inh_syn.x_S = 1
# ##############################################################################
# # Monitors
# ##############################################################################
# Note that we could use a single monitor for all neurons instead, but in this
# way plotting is a bit easier in the end
exc_mon = SpikeMonitor(exc_neurons)
inh_mon = SpikeMonitor(inh_neurons)
### We record some additional data from a single excitatory neuron
ni = 50
# Record conductances and membrane potential of neuron ni
state_mon = StateMonitor(exc_neurons, ['v', 'g_e', 'g_i'], record=ni)
# We make sure to monitor synaptic variables after synapse are updated in order
# to use simple recurrence relations to reconstruct them. Record all synapses
# originating from neuron ni
synapse_mon = StateMonitor(exc_syn, ['u_S', 'x_S'],
record=exc_syn[ni, :], when='after_synapses')
# ##############################################################################
# # Simulation run
# ##############################################################################
run(duration, report='text')
################################################################################
# Analysis and plotting
################################################################################
plt.style.use('figures.mplstyle')
### Spiking activity (w/ rate)
fig1, ax = plt.subplots(nrows=2, ncols=1, sharex=False,
gridspec_kw={'height_ratios': [3, 1],
'left': 0.18, 'bottom': 0.18, 'top': 0.95,
'hspace': 0.1},
figsize=(3.07, 3.07))
ax[0].plot(exc_mon.t[exc_mon.i <= N_e//4]/ms,
exc_mon.i[exc_mon.i <= N_e//4], '|', color='C0')
ax[0].plot(inh_mon.t[inh_mon.i <= N_i//4]/ms,
inh_mon.i[inh_mon.i <= N_i//4]+N_e//4, '|', color='C1')
pu.adjust_spines(ax[0], ['left'])
ax[0].set(xlim=(0., duration/ms), ylim=(0, (N_e+N_i)//4), ylabel='neuron index')
# Generate frequencies
bin_size = 1*ms
spk_count, bin_edges = np.histogram(np.r_[exc_mon.t/ms, inh_mon.t/ms],
int(duration/ms))
rate = double(spk_count)/(N_e + N_i)/bin_size/Hz
ax[1].plot(bin_edges[:-1], rate, '-', color='k')
pu.adjust_spines(ax[1], ['left', 'bottom'])
ax[1].set(xlim=(0., duration/ms), ylim=(0, 10.),
xlabel='time (ms)', ylabel='rate (Hz)')
pu.adjust_ylabels(ax, x_offset=-0.18)
### Dynamics of a single neuron
fig2, ax = plt.subplots(4, sharex=False,
gridspec_kw={'left': 0.27, 'bottom': 0.18, 'top': 0.95,
'hspace': 0.2},
figsize=(3.07, 3.07))
### Postsynaptic conductances
ax[0].plot(state_mon.t/ms, state_mon.g_e[0]/nS, color='C0')
ax[0].plot(state_mon.t/ms, -state_mon.g_i[0]/nS, color='C1')
ax[0].plot([state_mon.t[0]/ms, state_mon.t[-1]/ms], [0, 0], color='grey',
linestyle=':')
# Adjust axis
pu.adjust_spines(ax[0], ['left'])
ax[0].set(xlim=(0., duration/ms), ylim=(-5.0, 0.25),
ylabel=f"postsyn.\nconduct.\n(${sympy.latex(nS)}$)")
### Membrane potential
ax[1].axhline(V_th/mV, color='C2', linestyle=':') # Threshold
# Artificially insert spikes
ax[1].plot(state_mon.t/ms, state_mon.v[0]/mV, color='black')
ax[1].vlines(exc_mon.t[exc_mon.i == ni]/ms, V_th/mV, 0, color='black')
pu.adjust_spines(ax[1], ['left'])
ax[1].set(xlim=(0., duration/ms), ylim=(-1+V_r/mV, 0.),
ylabel=f"membrane\npotential\n(${sympy.latex(mV)}$)")
### Synaptic variables
# Retrieves indexes of spikes in the synaptic monitor using the fact that we
# are sampling spikes and synaptic variables by the same dt
spk_index = np.in1d(synapse_mon.t, exc_mon.t[exc_mon.i == ni])
ax[2].plot(synapse_mon.t[spk_index]/ms, synapse_mon.x_S[0][spk_index], '.',
ms=4, color='C3')
ax[2].plot(synapse_mon.t[spk_index]/ms, synapse_mon.u_S[0][spk_index], '.',
ms=4, color='C4')
# Super-impose reconstructed solutions
time = synapse_mon.t # time vector
tspk = Quantity(synapse_mon.t, copy=True) # Spike times
for ts in exc_mon.t[exc_mon.i == ni]:
tspk[time >= ts] = ts
ax[2].plot(synapse_mon.t/ms, 1 + (synapse_mon.x_S[0]-1)*exp(-(time-tspk)*Omega_d),
'-', color='C3')
ax[2].plot(synapse_mon.t/ms, synapse_mon.u_S[0]*exp(-(time-tspk)*Omega_f),
'-', color='C4')
# Adjust axis
pu.adjust_spines(ax[2], ['left'])
ax[2].set(xlim=(0., duration/ms), ylim=(-0.05, 1.05),
ylabel='synaptic\nvariables\n$u_S,\,x_S$')
nspikes = np.sum(spk_index)
x_S_spike = synapse_mon.x_S[0][spk_index]
u_S_spike = synapse_mon.u_S[0][spk_index]
ax[3].vlines(synapse_mon.t[spk_index]/ms, np.zeros(nspikes),
x_S_spike*u_S_spike/(1-u_S_spike))
pu.adjust_spines(ax[3], ['left', 'bottom'])
ax[3].set(xlim=(0., duration/ms), ylim=(-0.01, 0.62),
yticks=np.arange(0, 0.62, 0.2), xlabel='time (ms)', ylabel='$r_S$')
pu.adjust_ylabels(ax, x_offset=-0.20)
plt.show()
|
