1/*
2 *  brunel-2000_newconnect.sli
3 *
4 *  This file is part of NEST.
5 *
6 *  Copyright (C) 2004 The NEST Initiative
7 *
8 *  NEST is free software: you can redistribute it and/or modify
9 *  it under the terms of the GNU General Public License as published by
10 *  the Free Software Foundation, either version 2 of the License, or
11 *  (at your option) any later version.
12 *
13 *  NEST is distributed in the hope that it will be useful,
14 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
15 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
16 *  GNU General Public License for more details.
17 *
18 *  You should have received a copy of the GNU General Public License
19 *  along with NEST.  If not, see <http://www.gnu.org/licenses/>.
20 *
21 */
22
23% autorun=true
24
25/*
26   Brunel Network
27
28   The SLI code in this file creates a sparsely coupled network of
29   excitatory and inhibitory neurons.  Connections within and across
30   both populations are created at random.  Both neuron populations
31   receive Poissonian background input.  The spike output of 50
32   neurons from each population are recorded.  Neurons are modeled
33   as leaky integrate-and-fire neurons with current-injecting synapses
34   (alpha functions).  The model is based on
35
36      Nicolas Brunel
37      Dynamics of sparsely connected networks of excitatory
38      and inhibitory spiking neurons
39      Journal of Computational Neuroscience, 2000, vol 8, pp 183-208.
40
41   There are two differences to Brunel's model: we use alpha
42   functions instead of delta for synaptic currents, and our neurons
43   reset to the resting potential (0 mv) instead of to half-way
44   between resting potential and threshold.
45
46   This example shows how to
47
48      - create populations
49      - instrument a network with injection and recording devices
50      - record data to files
51      - define own functions
52      - set parameters in a simple way
53      - communicate with the user in a simple way
54
55   This brunel-2000_newconnect.sli version shows how to construct
56   the network using the new Connect function introduced with
57   NEST 2.4.
58
59   Abigail Morrison, Marc-Oliver Gewaltig, Hans Ekkehard Plesser
60*/
61
62%%% PARAMETER SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63
64% define all relevant parameters: changes should be made here
65% all data is place in the userdict dictionary
66
67/order 2500 def     % scales size of network (total 5*order neurons)
68
69% case C : asynchronous irregullar
70/g      5.0 def    % rel strength, inhibitory synapses
71/eta    2.0 def    % nu_ext / nu_thresh
72
73% case D : slow oscillations
74%/g      4.5  def    % rel strength, inhibitory synapses
75%/eta    0.95 def    % nu_ext / nu_thresh
76
77/simtime  1000.0 def % simulation time [ms]
78/dt          0.1 def % simulation step length [ms]
79
80
81%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
82%       NO USER-SERVICABLE PARTS BELOW
83%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84
85%%% FUNCTION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
86
87% Take spike recorder, find total number of spikes registered,
88% return average rate per neuron in Hz.
89% NOTE: If you are running with several MPI processes, this
90%       function only gives an approximation to the true rate.
91%
92% spike_det ComputeRate -> rate
93/ComputeRate
94{
95  << >> begin  % anonymous dictionary for local variables
96
97    /sr Set
98
99    % We need to guess how many neurons we record from.
100    % This assumes an even distribution of nodes across
101    % processes, as well as homogeneous activity in the
102    % network. So this is really a hack. NEST needs better
103    % support for rate calculations, such as giving the
104    % number of neurons recorded from by each spike recorder.
105
106    userdict /Nrec get cvd NumProcesses div /nnrn Set
107    sr /n_events get nnrn userdict /simtime get mul div
108    1000 mul         % convert from mHz to Hz, leave on stack
109
110  end
111} bind             % optional, improves performance
112def
113
114
115% Compute the maximum of postsynaptic potential
116% for a synaptic input current of unit amplitude
117% (1 pA)
118/ComputePSPnorm
119{
120  % calculate the normalization factor for the PSP
121  (
122  a = tauMem / tauSyn;
123  b = 1.0 / tauSyn - 1.0 / tauMem;
124  % time of maximum
125  t_max = 1.0/b * (-LambertWm1(-exp(-1.0/a)/a)-1.0/a);
126  % maximum of PSP for current of unit amplitude
127  exp(1.0)/(tauSyn*CMem*b) * ((exp(-t_max/tauMem) - exp(-t_max/tauSyn)) / b - t_max*exp(-t_max/tauSyn))
128  ) ExecMath
129}
130def
131
132%%% PREPARATION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
133
134/NE 4 order mul cvi def  % number of excitatory neurons
135/NI 1 order mul cvi def  % number of inhibitory neurons
136/N  NI NE add def        % total number of neurons
137
138/epsilon 0.1 def            % connectivity
139/CE epsilon NE mul cvi def  % number of excitatory synapses on neuron
140/CI epsilon NI mul cvi def  % number of inhibitory synapses on neuron
141/C  CE CI add def           % total number of internal synapses per n.
142/Cext CE def                % number of external synapses on neuron
143
144/tauMem 20.0 def    % neuron membrane time constant [ms]
145/CMem    1.0 def    % membrane capacity [pF]
146/tauSyn  0.5 def    % synaptic time constant [ms]
147/tauRef  2.0 def    % refractory time [ms]
148/E_L     0.0 def    % resting potential [mV]
149/theta  20.0 def    % threshold
150
151
152% amplitude of PSP given 1pA current
153ComputePSPnorm /J_max_unit Set
154
155% synaptic weights, scaled for our alpha functions, such that
156% for constant membrane potential, the peak amplitude of the PSP
157% equals J
158
159/delay   1.5 def         % synaptic delay, all connections [ms]
160/J       0.1 def         % synaptic weight [mV]
161/JE J J_max_unit div def % synaptic weight [pA]
162/JI g JE mul neg def     % inhibitory
163
164% threshold rate, equivalent rate of events needed to
165% have mean input current equal to threshold
166/nu_thresh ((theta * CMem) / (JE*CE*exp(1)*tauMem*tauSyn)) ExecMath def
167/nu_ext eta nu_thresh mul def     % external rate per synapse
168/p_rate nu_ext Cext mul 1000. mul def % external input rate per neuron
169                                        % must be given in Hz
170
171% number of neurons per population to record from
172/Nrec 50 def
173
174
175%%% CONSTRUCTION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
176
177ResetKernel      % clear all existing network elements
178M_WARNING setverbosity
179
180% Number of threads per program instance. When using MPI, the mpirun
181% call determines the number of MPI processes (= program instances).
182% The total number of virtual processes is
183%     #MPI processes x local_num_threads.
184%
185% The number of local_num_threads can be given to the script using
186% the --userargs commandline switch like this:
187%     nest --userargs=threads=4 script.sli
188%
189% If it is not given, it defaults to 2
190statusdict/userargs :: size 1 geq
191{
192  0 get (=) breakup 1 get int /local_num_threads Set
193}
194{
195  2 /local_num_threads Set
196} ifelse
197
198% set resolution and total/local number of threads
199<<
200    /resolution  dt
201    /local_num_threads local_num_threads
202    /overwrite_files true
203>> SetKernelStatus
204
205 tic % start timer on construction
206
207% Setting neuron default parameters
208
209      (Configuring neuron parameters.) =
210      /iaf_psc_alpha
211        <<
212          /tau_m       tauMem
213          /tau_syn_ex  tauSyn
214          /tau_syn_in  tauSyn
215          /t_ref       tauRef
216          /E_L     E_L
217          /V_m     E_L
218          /V_th    theta
219          /V_reset    10.0 % according to Brunel (2000) p 185
220          /C_m     1.0     % Yes, capacitance is 1 in the Brunel model
221        >> SetDefaults
222
223      (Creating the network.) =  % show message
224
225      /E_neurons /iaf_psc_alpha NE Create def
226      /I_neurons /iaf_psc_alpha NI Create def
227
228      % list of all neurons
229      /allNeurons E_neurons I_neurons join def
230
231      (Creating Poisson generators.) =
232      /poisson_generator
233      <<                % set firing rate
234        /rate p_rate
235      >> SetDefaults
236
237      /expoisson /poisson_generator Create def
238      /inpoisson /poisson_generator Create def
239
240      % We could do with only one recorder,
241      % but by recording the populations separately,
242      % we needn't sort the recordings later
243      (Creating excitatory spike recorder.) =
244      /exsr /spike_recorder Create def
245      exsr
246      <<
247        /label (brunel-2-ex-threaded)
248        /record_to /ascii
249      >> SetStatus
250
251      (Creating inhibitory spike recorder.) =
252      /insr /spike_recorder Create def
253      insr
254      <<
255        /label (brunel-2-in-threaded)
256        /record_to /ascii
257      >> SetStatus
258
259      % Create custom synapse types with appropriate values for
260      % our excitatory and inhibitory connections
261      /static_synapse_hom_w << /delay delay >> SetDefaults
262      /static_synapse_hom_w /syn_ex << /weight JE >> CopyModel
263      /static_synapse_hom_w /syn_in << /weight JI >> CopyModel
264
265      (Connecting excitatory neurons.) =
266      expoisson E_neurons /all_to_all /syn_ex Connect
267
268      E_neurons E_neurons << /rule /fixed_indegree /indegree CE >> << /synapse_model /syn_ex >> Connect
269
270      I_neurons E_neurons << /rule /fixed_indegree /indegree CI >> << /synapse_model /syn_in >> Connect
271
272      (Connecting inhibitory population.) =
273      inpoisson I_neurons /all_to_all /syn_ex Connect
274
275      E_neurons I_neurons << /rule /fixed_indegree /indegree CE >> << /synapse_model /syn_ex >> Connect
276
277      I_neurons I_neurons << /rule /fixed_indegree /indegree CI >> << /synapse_model /syn_in >> Connect
278
279      % Spike recorders are connected to the first Nrec neurons in each
280      % population.  Since neurons are equal and connectivity is homogeneously
281      % randomized, this is equivalent to picking Nrec neurons at random
282      % from each population
283      (Connecting spike recorders.) =
284
285      E_neurons Nrec Take exsr Connect  % pick the first 500 neurons
286      I_neurons Nrec Take insr Connect
287
288      toc /BuildCPUTime Set
289
290      /syn_ex GetDefaults /num_connections get /n_ex_syn Set
291      /syn_in GetDefaults /num_connections get /n_in_syn Set
292
293%%% SIMULATION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
294
295  % run, measure computer time with tic-toc
296  (Starting simulation.) =
297  tic
298  simtime Simulate
299  toc /SimCPUTime Set
300
301  % write a little report
302  (\nBrunel Network Simulation) =
303  (Number of Threads : ) =only local_num_threads =
304  (Number of Neurons : ) =only N =
305  (Number of Synapses: ) =only GetKernelStatus /num_connections get =
306  (       Excitatory : ) =only n_ex_syn =
307  (       Inhibitory : ) =only n_in_syn =
308  (Excitatory rate   : ) =only exsr ComputeRate =only ( Hz) =
309  (Inhibitory rate   : ) =only insr ComputeRate =only ( Hz) =
310  (Building time     : ) =only BuildCPUTime =only ( s) =
311  (Simulation time   : ) =only SimCPUTime   =only ( s\n) =
312
313