1# -*- coding: utf-8 -*-
2#
3# brunel_exp_multisynapse_nest.py
4#
5# This file is part of NEST.
6#
7# Copyright (C) 2004 The NEST Initiative
8#
9# NEST is free software: you can redistribute it and/or modify
10# it under the terms of the GNU General Public License as published by
11# the Free Software Foundation, either version 2 of the License, or
12# (at your option) any later version.
13#
14# NEST is distributed in the hope that it will be useful,
15# but WITHOUT ANY WARRANTY; without even the implied warranty of
16# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
17# GNU General Public License for more details.
18#
19# You should have received a copy of the GNU General Public License
20# along with NEST.  If not, see <http://www.gnu.org/licenses/>.
21
22"""
23Random balanced network (exp synapses, multiple time constants)
24---------------------------------------------------------------
25
26This script simulates an excitatory and an inhibitory population on
27the basis of the network used in [1]_.
28
29The example demonstrate the usage of the multisynapse neuron
30model. Each spike arriving at the neuron triggers an exponential
31PSP. The time constant associated with the PSP is defined in the
32receptor type array tau_syn of each neuron. The receptor types of all
33connections are uniformally distributed, resulting in uniformally
34distributed time constants of the PSPs.
35
36When connecting the network, customary synapse models are used, which
37allow for querying the number of created synapses. Using spike
38recorders, the average firing rates of the neurons in the populations
39are established. The building as well as the simulation time of the
40network are recorded.
41
42References
43~~~~~~~~~~
44
45.. [1] Brunel N (2000). Dynamics of sparsely connected networks of excitatory and
46       inhibitory spiking neurons. Journal of Computational Neuroscience 8,
47       183-208.
48
49See Also
50~~~~~~~~
51
52:doc:`brunel_alpha_nest`
53
54"""
55
56###############################################################################
57# Import all necessary modules for simulation, analysis and plotting.
58
59import time
60import nest
61import nest.raster_plot
62import matplotlib.pyplot as plt
63
64nest.ResetKernel()
65
66###############################################################################
67# Assigning the current time to a variable in order to determine the build
68# time of the network.
69
70startbuild = time.time()
71
72###############################################################################
73# Assigning the simulation parameters to variables.
74
75dt = 0.1    # the resolution in ms
76simtime = 1000.0  # Simulation time in ms
77delay = 1.5    # synaptic delay in ms
78
79###############################################################################
80# Definition of the parameters crucial for asynchronous irregular firing of
81# the neurons.
82
83g = 5.0  # ratio inhibitory weight/excitatory weight
84eta = 2.0  # external rate relative to threshold rate
85epsilon = 0.1  # connection probability
86
87###############################################################################
88# Definition of the number of neurons in the network and the number of neurons
89# recorded from
90
91order = 2500
92NE = 4 * order  # number of excitatory neurons
93NI = 1 * order  # number of inhibitory neurons
94N_neurons = NE + NI   # number of neurons in total
95N_rec = 50      # record from 50 neurons
96
97###############################################################################
98# Definition of connectivity parameters
99
100CE = int(epsilon * NE)  # number of excitatory synapses per neuron
101CI = int(epsilon * NI)  # number of inhibitory synapses per neuron
102C_tot = int(CI + CE)      # total number of synapses per neuron
103
104###############################################################################
105# Initialization of the parameters of the integrate and fire neuron and the
106# synapses. The parameters of the neuron are stored in a dictionary.
107
108tauMem = 20.0  # time constant of membrane potential in ms
109theta = 20.0  # membrane threshold potential in mV
110J = 0.1   # postsynaptic amplitude in mV
111nr_ports = 100  # number of receptor types
112# Create array of synaptic time constants for each neuron,
113# ranging from 0.1 to 1.09 ms.
114tau_syn = [0.1 + 0.01 * i for i in range(nr_ports)]
115neuron_params = {"C_m": 1.0,
116                 "tau_m": tauMem,
117                 "t_ref": 2.0,
118                 "E_L": 0.0,
119                 "V_reset": 0.0,
120                 "V_m": 0.0,
121                 "V_th": theta,
122                 "tau_syn": tau_syn}
123J_ex = J       # amplitude of excitatory postsynaptic current
124J_in = -g * J_ex  # amplitude of inhibitory postsynaptic current
125
126###############################################################################
127# Definition of threshold rate, which is the external rate needed to fix the
128# membrane potential around its threshold, the external firing rate and the
129# rate of the poisson generator which is multiplied by the in-degree CE and
130# converted to Hz by multiplication by 1000.
131
132nu_th = theta / (J * CE * tauMem)
133nu_ex = eta * nu_th
134p_rate = 1000.0 * nu_ex * CE
135
136###############################################################################
137# Configuration of the simulation kernel by the previously defined time
138# resolution used in the simulation. Setting ``print_time`` to `True` prints the
139# already processed simulation time as well as its percentage of the total
140# simulation time.
141
142nest.resolution = dt
143nest.print_time = True
144nest.overwrite_files = True
145nest.local_num_threads = 4
146
147print("Building network")
148
149###############################################################################
150# Creation of the nodes using ``Create``. We store the returned handles in
151# variables for later reference. Here the excitatory and inhibitory, as well
152# as the poisson generator and two spike recorders. The spike recorders will
153# later be used to record excitatory and inhibitory spikes. Properties of the
154# nodes are specified via ``params``, which expects a dictionary.
155
156nodes_ex = nest.Create("iaf_psc_exp_multisynapse", NE, params=neuron_params)
157nodes_in = nest.Create("iaf_psc_exp_multisynapse", NI, params=neuron_params)
158noise = nest.Create("poisson_generator", params={"rate": p_rate})
159espikes = nest.Create("spike_recorder")
160ispikes = nest.Create("spike_recorder")
161
162################################################################################
163# Configuration of the spike recorders recording excitatory and inhibitory
164# spikes by sending parameter dictionaries to ``set``. Setting the property
165# `record_to` to *"ascii"* ensures that the spikes will be recorded to a file,
166# whose name starts with the string assigned to the property `label`.
167
168espikes.set(label="brunel-py-ex", record_to="ascii")
169ispikes.set(label="brunel-py-in", record_to="ascii")
170
171print("Connecting devices")
172
173###############################################################################
174# Definition of a synapse using ``CopyModel``, which expects the model name of
175# a pre-defined synapse, the name of the customary synapse and an optional
176# parameter dictionary. The parameters defined in the dictionary will be the
177# default parameter for the customary synapse. Here we define one synapse for
178# the excitatory and one for the inhibitory connections giving the
179# previously defined weights and equal delays.
180
181nest.CopyModel("static_synapse", "excitatory",
182               {"weight": J_ex, "delay": delay})
183nest.CopyModel("static_synapse", "inhibitory",
184               {"weight": J_in, "delay": delay})
185
186###################################################################################
187# Connecting the previously defined poisson generator to the excitatory and
188# inhibitory neurons using the excitatory synapse. Since the poisson
189# generator is connected to all neurons in the population the default rule
190# (# ``all_to_all``) of ``Connect`` is used. The synaptic properties are
191# pre-defined # in a dictionary and inserted via ``syn_spec``. As synaptic model
192# the pre-defined synapses "excitatory" and "inhibitory" are choosen,
193# thus setting ``weight`` and ``delay``. The receptor type is drawn from a
194# distribution for each connection, which is specified in the synapse
195# properties by assigning a dictionary to the keyword ``receptor_type``,
196# which includes the specification of the distribution and the associated
197# parameter.
198
199syn_params_ex = {"synapse_model": "excitatory",
200                 "receptor_type": nest.random.uniform_int(max=nr_ports - 1) + 1
201                 }
202syn_params_in = {"synapse_model": "inhibitory",
203                 "receptor_type": nest.random.uniform_int(max=nr_ports - 1) + 1
204                 }
205
206nest.Connect(noise, nodes_ex, syn_spec=syn_params_ex)
207nest.Connect(noise, nodes_in, syn_spec=syn_params_ex)
208
209###############################################################################
210# Connecting the first ``N_rec`` nodes of the excitatory and inhibitory
211# population to the associated spike recorders using excitatory synapses.
212# Here the same shortcut for the specification of the synapse as defined
213# above is used.
214
215nest.Connect(nodes_ex[:N_rec], espikes, syn_spec="excitatory")
216nest.Connect(nodes_in[:N_rec], ispikes, syn_spec="excitatory")
217
218print("Connecting network")
219
220print("Excitatory connections")
221
222###############################################################################
223# Connecting the excitatory population to all neurons while distribution the
224# ports. Here we use the previously defined parameter dictionary
225# ``syn_params_ex``. Beforehand, the connection parameter are defined in a
226# dictionary. Here we use the connection rule ``fixed_indegree``,
227# which requires the definition of the indegree.
228
229conn_params_ex = {'rule': 'fixed_indegree', 'indegree': CE}
230nest.Connect(nodes_ex, nodes_ex + nodes_in, conn_params_ex, syn_params_ex)
231
232print("Inhibitory connections")
233
234###############################################################################
235# Connecting the inhibitory population to all neurons while distribution the
236# ports. Here we use the previously defined parameter dictionary
237# ``syn_params_in``.The connection parameter are defined analogously to the
238# connection from the excitatory population defined above.
239
240conn_params_in = {'rule': 'fixed_indegree', 'indegree': CI}
241nest.Connect(nodes_in, nodes_ex + nodes_in, conn_params_in, syn_params_in)
242
243###############################################################################
244# Storage of the time point after the buildup of the network in a variable.
245
246endbuild = time.time()
247
248###############################################################################
249# Simulation of the network.
250
251print("Simulating")
252
253nest.Simulate(simtime)
254
255###############################################################################
256# Storage of the time point after the simulation of the network in a variable.
257
258endsimulate = time.time()
259
260###############################################################################
261# Reading out the total number of spikes received from the spike recorder
262# connected to the excitatory population and the inhibitory population.
263
264events_ex = espikes.n_events
265events_in = ispikes.n_events
266
267###############################################################################
268# Calculation of the average firing rate of the excitatory and the inhibitory
269# neurons by dividing the total number of recorded spikes by the number of
270# neurons recorded from and the simulation time. The multiplication by 1000.0
271# converts the unit 1/ms to 1/s=Hz.
272
273rate_ex = events_ex / simtime * 1000.0 / N_rec
274rate_in = events_in / simtime * 1000.0 / N_rec
275
276###############################################################################
277# Reading out the number of connections established using the excitatory and
278# inhibitory synapse model. The numbers are summed up resulting in the total
279# number of synapses.
280
281num_synapses = (nest.GetDefaults("excitatory")["num_connections"]
282                + nest.GetDefaults("inhibitory")["num_connections"])
283
284###############################################################################
285# Establishing the time it took to build and simulate the network by taking
286# the difference of the pre-defined time variables.
287
288build_time = endbuild - startbuild
289sim_time = endsimulate - endbuild
290
291###############################################################################
292# Printing the network properties, firing rates and building times.
293
294print("Brunel network simulation (Python)")
295print(f"Number of neurons : {N_neurons}")
296print(f"Number of synapses: {num_synapses}")
297print(f"       Exitatory  : {int(CE * N_neurons) + N_neurons}")
298print(f"       Inhibitory : {int(CI * N_neurons)}")
299print(f"Excitatory rate   : {rate_ex:.2f} Hz")
300print(f"Inhibitory rate   : {rate_in:.2f} Hz")
301print(f"Building time     : {build_time:.2f} s")
302print(f"Simulation time   : {sim_time:.2f} s")
303
304###############################################################################
305# Plot a raster of the excitatory neurons and a histogram.
306
307nest.raster_plot.from_device(espikes, hist=True)
308plt.show()
309