/*
* brunel-sli_connect.sli
*
* This file is part of NEST.
*
* Copyright (C) 2004 The NEST Initiative
*
* NEST is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* NEST is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with NEST. If not, see <http://www.gnu.org/licenses/>.
*
*/
/*
Brunel Network
The SLI code in this file creates a sparsely coupled network of
excitatory and inhibitory neurons. Connections within and across
both populations are created at random. Both neuron populations
receive Poissonian background input. The spike output of 50
neurons from each population are recorded. Neurons are modeled
as leaky integrate-and-fire neurons with current-injecting synapses
(alpha functions). The model is based on
Nicolas Brunel
Dynamics of sparsely connected networks of excitatory
and inhibitory spiking neurons
Journal of Computational Neuroscience, 2000, vol 8, pp 183-208.
This example is for instruction purposes only.
For research, we recommend to use the PyNEST version of this model
in the pynest example folder, or brunel-2000.sli in this folder.
This example demonstrates how to
- create and use different synapse models
- construct a random network, using sli and the plain connect
command alone
- instrument a network with injection and recording devices
- record data to files
- define own functions
- set parameters in a simple way
- communicate with the user in a simple way
- inline constants and functions to increase execution speed
Connections are drawn by one of the functions
build_network_sli_1
build_network_sli_2
build_network_sli_3
With increasing number, the functions are more optimised compared to its
predecessor.
Version 1 uses straight forward SLI commands.
Version 2 uses a variant of connect that expects the index of the synapse
model rather than its name.
Version 3 uses inlining and early binding of names and symbols in the inner
loops.
The version is selected by the variable /build_network, defined below.
Its default is 3.
Marc-Oliver Gewaltig, based on brunel-2000.sli by
Abigail Morrison, Marc-Oliver Gewaltig and Hans Ekkehard Plesser
*/
%%% PARAMETER SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% define all relevant parameters: changes should be made here
% all data is place in the userdict dictionary
/order 2500 def % scales size of network (total 5*order neurons)
% case C : asynchronous irregullar
/g 5.0 def % rel strength, inhibitory synapses
/eta 2.0 def % nu_ext / nu_thresh
% case D : slow oscillations
%/g 4.5 def % rel strength, inhibitory synapses
%/eta 0.95 def % nu_ext / nu_thresh
/simtime 1000.0 def % simulation time [ms]
/dt 0.1 def % simulation step length [ms]
% Number of POSIX threads per program instance.
% When using MPI, the mpirun call determines the number
% of MPI processes (=program instances). The total number
% of virtual processes is #MPI processes x local_num_threads.
/local_num_threads 8 def
% Which connect routine to use
% 1 simplest and slowest, 3 the fastest and most optimised.
/build_network 3 def
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% NO USER-SERVICABLE PARTS BELOW
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%% FUNCTION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Take spike detector, find total number of spikes registered,
% return average rate per neuron in Hz.
% NOTE: If you are running with several MPI processes, this
% function only gives an approximation to the true rate.
%
% spike_det ComputeRate -> rate
/ComputeRate
{
<< >> begin % anonymous dictionary for local variables
/sdet Set
% We need to guess how many neurons we record from.
% This assumes an even distribution of nodes across
% processes, as well as homogeneous activity in the
% network. So this is really a hack. NEST needs better
% support for rate calculations, such as giving the
% number of neurons recorded from by each spike detector.
userdict /Nrec get cvd NumProcesses div /nnrn Set
sdet /n_events get nnrn userdict /simtime get mul div
1000 mul % convert from mHz to Hz, leave on stack
end
} bind % optional, improves performance
def
% Compute the maximum of postsynaptic potential
% for a synaptic input current of unit amplitude
% (1 pA)
/ComputePSPnorm
{
% calculate the normalization factor for the PSP
(
a = tauMem / tauSyn;
b = 1.0 / tauSyn - 1.0 / tauMem;
% time of maximum
t_max = 1.0/b * (-LambertWm1(-exp(-1.0/a)/a)-1.0/a);
% maximum of PSP for current of unit amplitude
exp(1.0)/(tauSyn*CMem*b) * ((exp(-t_max/tauMem) - exp(-t_max/tauSyn)) / b - t_max*exp(-t_max/tauSyn))
) ExecMath
}
def
%%% PREPARATION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
/NE 4 order mul cvi def % number of excitatory neurons
/NI 1 order mul cvi def % number of inhibitory neurons
/N NI NE add def % total number of neurons
/epsilon 0.1 def % connectivity
/CE epsilon NE mul cvi def % number of excitatory synapses on neuron
/CI epsilon NI mul cvi def % number of inhibitory synapses on neuron
/C CE CI add def % total number of internal synapses per n.
/Cext CE def % number of external synapses on neuron
/tauMem 20.0 def % neuron membrane time constant [ms]
/CMem 1.0 def % membrane capacity [pF]
/tauSyn 0.5 def % synaptic time constant [ms]
/tauRef 2.0 def % refractory time [ms]
/U0 0.0 def % resting potential [mV]
/theta 20.0 def % threshold
% amplitude of PSP given 1pA current
ComputePSPnorm /J_max_unit Set
% synaptic weights, scaled for our alpha functions, such that
% for constant membrane potential, the peak amplitude of the PSP
% equals J
/delay 1.5 def % synaptic delay, all connections [ms]
/J 0.1 def % synaptic weight [mV]
/JE J J_max_unit div def % synaptic weight [pA]
/JI g JE mul neg def % inhibitory
% threshold rate, equivalent rate of events needed to
% have mean input current equal to threshold
/nu_thresh ((theta * CMem) / (JE*CE*exp(1)*tauMem*tauSyn)) ExecMath def
/nu_ext eta nu_thresh mul def % external rate per synapse
/p_rate nu_ext Cext mul 1000. mul def % external input rate per neuron
% must be given in Hz
% number of neurons per population to record from
/Nrec 50 def
/model /iaf_psc_alpha def
%%% CONSTRUCTION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
ResetKernel % clear all existing network elements
M_WARNING setverbosity
% set resolution and total/local number of threads
[0]
<<
/resolution dt
/local_num_threads local_num_threads
/print_time true
/overwrite_files true
>> SetStatus
% Setting neuron default parameters
(Configuring neuron parameters.) =
/iaf_psc_alpha
<<
/tau_m tauMem
/tau_syn_ex tauSyn
/tau_syn_in tauSyn
/t_ref tauRef
/E_L U0
/V_m U0
/V_th theta
/V_reset 10.0 % according to Brunel (2000) p 185
/C_m 1.0 % Yes, capacitance is 1 in the Brunel model
>> SetDefaults
(SLI only network building)
% Setting neuron default parameters
(Creating the network.) =
model NE Create /E_last Set
model NI Create /I_last Set
E_last 1 add /I_first Set
cout (First E neuron: ) <- 1 <- endl
(Last E neuron: ) <- E_last <- endl ;
cout (First I neuron: ) <- I_first <- endl
(Last I neuron: ) <- I_last <- endl ;
(Creating Poisson generators.) =
/poisson_generator
<< % set firing rate
/rate p_rate
>> SetDefaults
/expoisson /poisson_generator Create def
/inpoisson /poisson_generator Create def
% one detector would in principle be enough,
% but by recording the populations separately,
% we save us a lot of sorting work later
(Creating excitatory spike detector.) =
/exsd /spike_detector Create def
exsd
<<
/label (brunel-1-ex)
/withtime true % record time of spikes
/withgid true % record which neuron spiked
/to_file true % write results to a file
/to_memory false
>> SetStatus
(Creating inhibitory spike detector.) =
/insd /spike_detector Create def
insd
<<
/label (brunel-1-in)
/withtime true
/withgid true
/to_file true
/to_memory false
>> SetStatus
% Create custom synapse types with appropriate values for
% our excitatory and inhibitory connections
/static_synapse_hom_wd << /delay delay >> SetDefaults
/static_synapse_hom_wd /syn_ex << /weight JE >> CopyModel
/static_synapse_hom_wd /syn_in << /weight JI >> CopyModel
synapsedict /syn_ex get /syn_ex_id Set
synapsedict /syn_in get /syn_in_id Set
synapsedict /static_synapse get /syn_default Set
rngdict /MT19937 get 123456789 CreateRNG /rng Set
% This is the straight forward way to connect the Brunel network.
/build_network_sli_1
{
(Connecting background.) =
% Connect Poisson generator to neurons
1 1 E_last
{
expoisson exch /syn_ex Connect
} for
% Connect Poisson generator to neurons
I_first 1 I_last
{
inpoisson exch /syn_ex Connect
} for
% Spike detectors are connected to the first Nrec neurons in each
% population. Since neurons are equal and connectivity is homogeneously
% randomized, this is equivalent to picking Nrec neurons at random
% from each population
(Connecting spike detectors.) =
1 1 Nrec
{
exsd /syn_ex Connect
} for
I_first dup 1 exch Nrec add
{
insd /syn_ex Connect
} for
(Drawing excitatory connections.) =
% random convergent connect to all neurons
1 1 E_last
{
/target Set
% We have CE E -> E connections
CE
{
rng NE irand % draw random number in [0,NE)
1 add % plus 1 to get number in [1,NE]
target
/syn_ex Connect
} repeat
} for
I_first 1 I_last
{
/target Set
% E -> I connections
CE
{
rng NE irand 1 add
target /syn_ex Connect
} repeat
} for
(Drawing inhibitory connections) =
1 1 E_last
{
/target Set
% We have CI I -> E connections
CI
{
rng NI irand I_first add
target /syn_in Connect
} repeat
} for
% We must undef target to prevent the old value from
% being inlined
userdict /target undef
I_first 1 I_last
{
/target Set
CI
{
rng NI irand I_first add
target /syn_in Connect
} repeat
} for
} bind def
/build_network_sli_2
{
(Connecting background.) =
% Connect Poisson generator to neurons
1 1 E_last
{
expoisson exch syn_ex_id Connect_i_i_i
} for
% Connect Poisson generator to neurons
I_first 1 I_last
{
inpoisson exch syn_ex_id Connect_i_i_i
} for
% Spike detectors are connected to the first Nrec neurons in each
% population. Since neurons are equal and connectivity is homogeneously
% randomized, this is equivalent to picking Nrec neurons at random
% from each population
(Connecting spike detectors.) =
1 1 Nrec
{
exsd syn_ex_id Connect_i_i_i
} for
I_first dup 1 exch Nrec add
{
insd syn_ex_id Connect_i_i_i
} for
(Drawing excitatory connections.) =
% random convergent connect to all neurons
1 1 E_last
{
/target Set
% We have CE E -> E connections
CE
{
rng NE irand % draw random number in [0,NE)
1 add % plus 1 to get number in [1,NE]
target
syn_ex_id Connect_i_i_i
} repeat
} for
I_first 1 I_last
{
/target Set
% E -> I connections
CE
{
rng NE irand 1 add
target syn_ex_id Connect_i_i_i
} repeat
} for
(Drawing inhibitory connections) =
1 1 E_last
{
/target Set
% We have CI I -> E connections
CI
{
rng NI irand I_first add
target syn_in_id Connect_i_i_i
} repeat
} for
% We must undef target to prevent the old value from
% being inlined
userdict /target undef
I_first 1 I_last
{
/target Set
CI
{
rng NI irand I_first add
target syn_in_id Connect_i_i_i
} repeat
} for
} bind def
% This version squeezes the last bit of performance out of SLI
% - we inline all constants
% - we bind all command names
% The performance gain of this optimizations is around 20 %
/build_network_sli_3
{
(Connecting background.) =
% Connect Poisson generator to neurons
1 1 E_last
{
expoisson exch syn_ex_id Connect_i_i_i
} for
% Connect Poisson generator to neurons
I_first 1 I_last
{
inpoisson exch syn_ex_id Connect_i_i_i
} for
% Spike detectors are connected to the first Nrec neurons in each
% population. Since neurons are equal and connectivity is homogeneously
% randomized, this is equivalent to picking Nrec neurons at random
% from each population
(Connecting spike detectors.) =
1 1 Nrec
{
exsd syn_ex_id Connect_i_i_i
} for
I_first dup 1 exch Nrec add
{
insd syn_ex_id Connect_i_i_i
} for
(Drawing excitatory connections.) =
% random convergent connect to all neurons
1 1 E_last
{
/target Set_
% We have CE E -> E connections
CE
{
rng NE irand_g_i % draw random number in [0,NE)
1 add_ii % plus 1 to get number in [1,NE]
target
syn_ex_id Connect_i_i_i
} repeat_
} for
% We must undef target to prevent the old value from
% being inlined
userdict /target undef
I_first 1 I_last
{
/target Set_
% E -> I connections
CE
{
rng NE irand_g_i 1 add_ii
target syn_ex_id Connect_i_i_i
} repeat_
} for
(Drawing inhibitory connections) =
% We must undef target to prevent the old value from
% being inlined
userdict /target undef
1 1 E_last
{
/target Set_
% We have CI I -> E connections
CI
{
rng NI irand_g_i I_first add_ii
target syn_in_id Connect_i_i_i
} repeat_
} for
% We must undef target to prevent the old value from
% being inlined
userdict /target undef
I_first 1 I_last
{
/target Set_
CI
{
rng NI irand_g_i I_first add_ii
target syn_in_id Connect_i_i_i
} repeat_
} for
} bind userdict Inline def
%%% SIMULATION SECTION %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
tic
[
/build_network_sli_1
/build_network_sli_2
/build_network_sli_3
]
build_network 1 sub get
dup
(\n Builting with: ) =only ==
load exec
toc /BuildCPUTime Set
(Building time : ) =only BuildCPUTime =only ( s) =
/syn_ex GetDefaults /num_connections get /n_ex_syn Set
/syn_in GetDefaults /num_connections get /n_in_syn Set
% run, measure computer time with tic-toc
tic
simtime Simulate
toc /SimCPUTime Set
% write a little report
(\nBrunel Network Simulation) =
(Number of Neurons : ) =only N =
(Number of Synapses: ) =only 0 GetStatus /num_connections get =
( Excitatory : ) =only n_ex_syn =
( Inhibitory : ) =only n_in_syn =
(Excitatory rate : ) =only exsd ComputeRate =only ( Hz) =
(Inhibitory rate : ) =only insd ComputeRate =only ( Hz) =
(Building time : ) =only BuildCPUTime =only ( s) =
(Simulation time : ) =only SimCPUTime =only ( s\n) =