1# -*- coding: utf-8 -*- 2# 3# brunel_alpha_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 (alpha synapses) connected with NEST 24------------------------------------------------------------ 25 26This script simulates an excitatory and an inhibitory population on 27the basis of the network used in [1]_. 28 29In contrast to ``brunel-alpha-numpy.py``, this variant uses NEST's builtin 30connection routines to draw the random connections instead of NumPy. 31 32When connecting the network, customary synapse models are used, which 33allow for querying the number of created synapses. Using spike 34recorders, the average firing rates of the neurons in the populations 35are established. The building as well as the simulation time of the 36network are recorded. 37 38References 39~~~~~~~~~~ 40 41.. [1] Brunel N (2000). Dynamics of sparsely connected networks of excitatory and 42 inhibitory spiking neurons. Journal of Computational Neuroscience 8, 43 183-208. 44 45See Also 46~~~~~~~~ 47 48:doc:`brunel_alpha_numpy` 49 50""" 51 52############################################################################### 53# Import all necessary modules for simulation, analysis and plotting. Scipy 54# should be imported before nest. 55 56import time 57import numpy as np 58import scipy.special as sp 59 60import nest 61import nest.raster_plot 62import matplotlib.pyplot as plt 63 64 65############################################################################### 66# Definition of functions used in this example. First, define the `Lambert W` 67# function implemented in SLI. The second function computes the maximum of 68# the postsynaptic potential for a synaptic input current of unit amplitude 69# (1 pA) using the `Lambert W` function. Thus function will later be used to 70# calibrate the synaptic weights. 71 72 73def LambertWm1(x): 74 # Using scipy to mimic the gsl_sf_lambert_Wm1 function. 75 return sp.lambertw(x, k=-1 if x < 0 else 0).real 76 77 78def ComputePSPnorm(tauMem, CMem, tauSyn): 79 a = (tauMem / tauSyn) 80 b = (1.0 / tauSyn - 1.0 / tauMem) 81 82 # time of maximum 83 t_max = 1.0 / b * (-LambertWm1(-np.exp(-1.0 / a) / a) - 1.0 / a) 84 85 # maximum of PSP for current of unit amplitude 86 return (np.exp(1.0) / (tauSyn * CMem * b) * 87 ((np.exp(-t_max / tauMem) - np.exp(-t_max / tauSyn)) / b - 88 t_max * np.exp(-t_max / tauSyn))) 89 90 91nest.ResetKernel() 92 93############################################################################### 94# Assigning the current time to a variable in order to determine the build 95# time of the network. 96 97startbuild = time.time() 98 99 100############################################################################### 101# Assigning the simulation parameters to variables. 102 103dt = 0.1 # the resolution in ms 104simtime = 1000.0 # Simulation time in ms 105delay = 1.5 # synaptic delay in ms 106 107############################################################################### 108# Definition of the parameters crucial for asynchronous irregular firing of 109# the neurons. 110 111g = 5.0 # ratio inhibitory weight/excitatory weight 112eta = 2.0 # external rate relative to threshold rate 113epsilon = 0.1 # connection probability 114 115############################################################################### 116# Definition of the number of neurons in the network and the number of neurons 117# recorded from 118 119order = 2500 120NE = 4 * order # number of excitatory neurons 121NI = 1 * order # number of inhibitory neurons 122N_neurons = NE + NI # number of neurons in total 123N_rec = 50 # record from 50 neurons 124 125############################################################################### 126# Definition of connectivity parameters 127 128CE = int(epsilon * NE) # number of excitatory synapses per neuron 129CI = int(epsilon * NI) # number of inhibitory synapses per neuron 130C_tot = int(CI + CE) # total number of synapses per neuron 131 132############################################################################### 133# Initialization of the parameters of the integrate and fire neuron and the 134# synapses. The parameters of the neuron are stored in a dictionary. The 135# synaptic currents are normalized such that the amplitude of the PSP is J. 136 137tauSyn = 0.5 # synaptic time constant in ms 138tauMem = 20.0 # time constant of membrane potential in ms 139CMem = 250.0 # capacitance of membrane in in pF 140theta = 20.0 # membrane threshold potential in mV 141neuron_params = {"C_m": CMem, 142 "tau_m": tauMem, 143 "tau_syn_ex": tauSyn, 144 "tau_syn_in": tauSyn, 145 "t_ref": 2.0, 146 "E_L": 0.0, 147 "V_reset": 0.0, 148 "V_m": 0.0, 149 "V_th": theta} 150J = 0.1 # postsynaptic amplitude in mV 151J_unit = ComputePSPnorm(tauMem, CMem, tauSyn) 152J_ex = J / J_unit # amplitude of excitatory postsynaptic current 153J_in = -g * J_ex # amplitude of inhibitory postsynaptic current 154 155############################################################################### 156# Definition of threshold rate, which is the external rate needed to fix the 157# membrane potential around its threshold, the external firing rate and the 158# rate of the poisson generator which is multiplied by the in-degree CE and 159# converted to Hz by multiplication by 1000. 160 161nu_th = (theta * CMem) / (J_ex * CE * np.exp(1) * tauMem * tauSyn) 162nu_ex = eta * nu_th 163p_rate = 1000.0 * nu_ex * CE 164 165################################################################################ 166# Configuration of the simulation kernel by the previously defined time 167# resolution used in the simulation. Setting ``print_time`` to `True` prints the 168# already processed simulation time as well as its percentage of the total 169# simulation time. 170 171nest.resolution = dt 172nest.print_time = True 173nest.overwrite_files = True 174 175print("Building network") 176 177############################################################################### 178# Creation of the nodes using ``Create``. We store the returned handles in 179# variables for later reference. Here the excitatory and inhibitory, as well 180# as the poisson generator and two spike recorders. The spike recorders will 181# later be used to record excitatory and inhibitory spikes. Properties of the 182# nodes are specified via ``params``, which expects a dictionary. 183 184nodes_ex = nest.Create("iaf_psc_alpha", NE, params=neuron_params) 185nodes_in = nest.Create("iaf_psc_alpha", NI, params=neuron_params) 186noise = nest.Create("poisson_generator", params={"rate": p_rate}) 187espikes = nest.Create("spike_recorder") 188ispikes = nest.Create("spike_recorder") 189 190############################################################################### 191# Configuration of the spike recorders recording excitatory and inhibitory 192# spikes by sending parameter dictionaries to ``set``. Setting the property 193# `record_to` to *"ascii"* ensures that the spikes will be recorded to a file, 194# whose name starts with the string assigned to the property `label`. 195 196espikes.set(label="brunel-py-ex", record_to="ascii") 197ispikes.set(label="brunel-py-in", record_to="ascii") 198 199print("Connecting devices") 200 201############################################################################### 202# Definition of a synapse using ``CopyModel``, which expects the model name of 203# a pre-defined synapse, the name of the customary synapse and an optional 204# parameter dictionary. The parameters defined in the dictionary will be the 205# default parameter for the customary synapse. Here we define one synapse for 206# the excitatory and one for the inhibitory connections giving the 207# previously defined weights and equal delays. 208 209nest.CopyModel("static_synapse", "excitatory", 210 {"weight": J_ex, "delay": delay}) 211nest.CopyModel("static_synapse", "inhibitory", 212 {"weight": J_in, "delay": delay}) 213 214################################################################################# 215# Connecting the previously defined poisson generator to the excitatory and 216# inhibitory neurons using the excitatory synapse. Since the poisson 217# generator is connected to all neurons in the population the default rule 218# (``all_to_all``) of ``Connect`` is used. The synaptic properties are inserted 219# via ``syn_spec`` which expects a dictionary when defining multiple variables or 220# a string when simply using a pre-defined synapse. 221 222nest.Connect(noise, nodes_ex, syn_spec="excitatory") 223nest.Connect(noise, nodes_in, syn_spec="excitatory") 224 225############################################################################### 226# Connecting the first ``N_rec`` nodes of the excitatory and inhibitory 227# population to the associated spike recorders using excitatory synapses. 228# Here the same shortcut for the specification of the synapse as defined 229# above is used. 230 231nest.Connect(nodes_ex[:N_rec], espikes, syn_spec="excitatory") 232nest.Connect(nodes_in[:N_rec], ispikes, syn_spec="excitatory") 233 234print("Connecting network") 235 236print("Excitatory connections") 237 238############################################################################### 239# Connecting the excitatory population to all neurons using the pre-defined 240# excitatory synapse. Beforehand, the connection parameter are defined in a 241# dictionary. Here we use the connection rule ``fixed_indegree``, 242# which requires the definition of the indegree. Since the synapse 243# specification is reduced to assigning the pre-defined excitatory synapse it 244# suffices to insert a string. 245 246conn_params_ex = {'rule': 'fixed_indegree', 'indegree': CE} 247nest.Connect(nodes_ex, nodes_ex + nodes_in, conn_params_ex, "excitatory") 248 249print("Inhibitory connections") 250 251############################################################################### 252# Connecting the inhibitory population to all neurons using the pre-defined 253# inhibitory synapse. The connection parameter as well as the synapse 254# parameter are defined analogously to the connection from the excitatory 255# population defined above. 256 257conn_params_in = {'rule': 'fixed_indegree', 'indegree': CI} 258nest.Connect(nodes_in, nodes_ex + nodes_in, conn_params_in, "inhibitory") 259 260############################################################################### 261# Storage of the time point after the buildup of the network in a variable. 262 263endbuild = time.time() 264 265############################################################################### 266# Simulation of the network. 267 268print("Simulating") 269 270nest.Simulate(simtime) 271 272############################################################################### 273# Storage of the time point after the simulation of the network in a variable. 274 275endsimulate = time.time() 276 277############################################################################### 278# Reading out the total number of spikes received from the spike recorder 279# connected to the excitatory population and the inhibitory population. 280 281events_ex = espikes.n_events 282events_in = ispikes.n_events 283 284############################################################################### 285# Calculation of the average firing rate of the excitatory and the inhibitory 286# neurons by dividing the total number of recorded spikes by the number of 287# neurons recorded from and the simulation time. The multiplication by 1000.0 288# converts the unit 1/ms to 1/s=Hz. 289 290rate_ex = events_ex / simtime * 1000.0 / N_rec 291rate_in = events_in / simtime * 1000.0 / N_rec 292 293############################################################################### 294# Reading out the number of connections established using the excitatory and 295# inhibitory synapse model. The numbers are summed up resulting in the total 296# number of synapses. 297 298num_synapses = (nest.GetDefaults("excitatory")["num_connections"] + 299 nest.GetDefaults("inhibitory")["num_connections"]) 300 301############################################################################### 302# Establishing the time it took to build and simulate the network by taking 303# the difference of the pre-defined time variables. 304 305build_time = endbuild - startbuild 306sim_time = endsimulate - endbuild 307 308############################################################################### 309# Printing the network properties, firing rates and building times. 310 311print("Brunel network simulation (Python)") 312print(f"Number of neurons : {N_neurons}") 313print(f"Number of synapses: {num_synapses}") 314print(f" Exitatory : {int(CE * N_neurons) + N_neurons}") 315print(f" Inhibitory : {int(CI * N_neurons)}") 316print(f"Excitatory rate : {rate_ex:.2f} Hz") 317print(f"Inhibitory rate : {rate_in:.2f} Hz") 318print(f"Building time : {build_time:.2f} s") 319print(f"Simulation time : {sim_time:.2f} s") 320 321############################################################################### 322# Plot a raster of the excitatory neurons and a histogram. 323 324nest.raster_plot.from_device(espikes, hist=True) 325plt.show() 326