A sparsely connected recurrent network using Nest¶
In the last case study,
we use Uncertainpy to perform a feature based analysis of the sparsely
connected recurrent network by Brunel (2000).
We implement the Brunel network using NEST inside a Python function,
and create \(10000\) inhibitory and \(2500\) excitatory neurons.
We record the output from \(20\) of the excitatory neurons,
and simulate the network for \(1000\) ms.
This is the values used to create the results in the Uncertainpy paper.
If you want to just test the network, we recommend reducing the model to
\(2000\) inhibitory and \(500\) excitatory neurons,
and only simulate the network for \(100\) ms.
To be able to run this example you require NEST to be anle to run the model and
elephant, neo, and quantities to be able to use the network features.
We want to use NestModel to create our model.
NestModel requires the model function to be specified through
the run argument, unlike NeuronModel.
The NEST model function has the same requirements as a regular model function,
except it is restricted to return only two objects:
the final simulation time (simulation_end),
and a list of spike times for each neuron in the network (spiketrains).
NestModel then postproccess this result for us to a regular result.
The final uncertainty quantification of a NEST network therefore predicts the
probability for a spike to occur at any specific time point in the simulation.
We implement the Brunel network as such a function
(found in /examples/brunel/brunel.py):
import nest
def brunel_network(eta=2, g=2, delay=1.5, J=0.1):
"""
A sparsely connected recurrent network (Brunel).
Brunel N, Dynamics of Sparsely Connected Networks of Excitatory and
Inhibitory Spiking Neurons, Journal of Computational Neuroscience 8,
183-208 (2000).
Parameters
----------
eta : {int, float}, optional
External rate relative to threshold rate. Default is 2.
g : {int, float}, optional
Ratio inhibitory weight/excitatory weight. Default is 5.
delay : {int, float}, optional
Synaptic delay in ms. Default is 1.5.
J : {int, float}, optional
Amplitude of excitatory postsynaptic current. Default is 0.1
Notes
-----
Brunel N, Dynamics of Sparsely Connected Networks of Excitatory and
Inhibitory Spiking Neurons, Journal of Computational Neuroscience 8,
183-208 (2000).
"""
# Network parameters
N_rec = 20 # Record from 20 neurons
simulation_end = 1000 # Simulation time
tau_m = 20.0 # Time constant of membrane potential in ms
V_th = 20.0
N_E = 10000 # Number of excitatory neurons
N_I = 2500 # Number of inhibitory neurons
N_neurons = N_E + N_I # Number of neurons in total
C_E = int(N_E/10) # Number of excitatory synapses per neuron
C_I = int(N_I/10) # Number of inhibitory synapses per neuron
J_I = -g*J # Amplitude of inhibitory postsynaptic current
cutoff = 100 # Cutoff to avoid transient effects, in ms
nu_ex = eta*V_th/(J*C_E*tau_m)
p_rate = 1000.0*nu_ex*C_E
nest.ResetKernel()
# Configure kernel
nest.SetKernelStatus({"grng_seed": 10})
nest.SetDefaults('iaf_psc_delta',
{'C_m': 1.0,
'tau_m': tau_m,
't_ref': 2.0,
'E_L': 0.0,
'V_th': V_th,
'V_reset': 10.0})
# Create neurons
nodes = nest.Create('iaf_psc_delta', N_neurons)
nodes_E = nodes[:N_E]
nodes_I = nodes[N_E:]
noise = nest.Create('poisson_generator',1,{'rate': p_rate})
spikes = nest.Create('spike_detector',2,
[{'label': 'brunel-py-ex'},
{'label': 'brunel-py-in'}])
spikes_E = spikes[:1]
spikes_I = spikes[1:]
# Connect neurons to each other
nest.CopyModel('static_synapse_hom_w', 'excitatory',
{'weight':J, 'delay':delay})
nest.Connect(nodes_E, nodes,
{'rule': 'fixed_indegree', 'indegree': C_E},
'excitatory')
nest.CopyModel('static_synapse_hom_w', 'inhibitory',
{'weight': J_I, 'delay': delay})
nest.Connect(nodes_I, nodes,
{'rule': 'fixed_indegree', 'indegree': C_I},
'inhibitory')
# Connect poisson generator to all nodes
nest.Connect(noise, nodes, syn_spec='excitatory')
nest.Connect(nodes_E[:N_rec], spikes_E)
nest.Connect(nodes_I[:N_rec], spikes_I)
# Run the simulation
nest.Simulate(simulation_end)
events_E = nest.GetStatus(spikes_E, 'events')[0]
events_I = nest.GetStatus(spikes_I, 'events')[0]
# Excitatory spike trains
# Makes sure the spiketrain is added even if there are no results
# to get a regular result
spiketrains = []
for sender in nodes_E[:N_rec]:
spiketrain = events_E["times"][events_E["senders"] == sender]
spiketrain = spiketrain[spiketrain > cutoff] - cutoff
spiketrains.append(spiketrain)
simulation_end -= cutoff
return simulation_end, spiketrains
And use it to create our model (example found in
/examples/brunel/uq_brunel.py):
We set ignore=True since we are not interested in
the model result itself.
This is recommended for NEST models as long as you do not
need the model results, since the uncertainty calculations for the
for the model results require much time and memory.
# Create a Nest model from the brunel network function
# We set ``ignore=True`` since we are not interested in
# the model result itself.
# This is recommended for NEST models as long as you do not
# need the model results, since the uncertainty calculations for the
# for the model results require much time and memory.
model = un.NestModel(run=brunel_network, ignore=True)
The Brunel model has four uncertain parameters:
- the external rate (\(\nu_\mathrm{ext}\)) relative to threshold rate (\(\nu_\mathrm{thr}\)) given as \(\eta = \nu_\mathrm{ext}/\nu_\mathrm{thr}\),
- the relative strength of the inhibitory synapses \(g\),
- the synaptic delay \(D\), and
- the amplitude of excitatory postsynaptic current \(J_e\).
Depending on the parameterizations of the model, the Brunel network may be in several different activity states. For the current example, we limit our analysis to two of these states. We create two sets of parameters, one for each of two states, and assume the parameter uncertainties are characterized by uniform probability distributions within the ranges below:
| Parameter | Range SR | Range AI | Variable | Meaning |
|---|---|---|---|---|
| \(\eta\) | \([1.5, 3.5]\) | \([1.5, 3.5]\) | eta |
External rate relative to threshold rate |
| \(g\) | \([1, 3]\) | \([5, 8]\) | g |
Relative strength of inhibitory synapses |
| \(D\) | \([1.5, 3]\) | \([1.5, 3]\) | delay |
Synaptic delay (ms) |
| \(J_e\) | \([0.05, 0.15]\) | \([0.05, 0.15]\) | J_e |
Amplitude excitatory postsynaptic current (mV) |
These ranges correspond to the synchronous regular (SR) state, where the neurons are almost completely synchronized, and the asynchronous irregular (AI) state, where the neurons fire individually at low rates. We create two sets of parameters, one for each state:
# Parametes for the synchronous regular (SR) state
parameters = {"eta": cp.Uniform(1.5, 3.5),
"g": cp.Uniform(1, 3),
"delay": cp.Uniform(1.5, 3)}
parameters_SR = un.Parameters(parameters)
# Parameter for the asynchronous irregular (AI) state
parameters = {"eta": cp.Uniform(1.5, 2.2),
"g": cp.Uniform(5, 8),
"delay": cp.Uniform(1.5, 3)}
parameters_AI = un.Parameters(parameters)
We use the features in NetworkFeatures to examine features of the Brunel network.
features = un.NetworkFeatures()
We set up the problems with the SR parameter set and use polynomial chaos with point collocation to perform the uncertainty quantification and sensitivity analysis. We specify a filename for the data, and a folder where to save the figures, to keep the results from the AI and SR state separated. We also set the seed to easier be able to reproduce the result.
UQ = un.UncertaintyQuantification(model,
parameters=parameters_SR,
features=features)
# Perform uncertainty quantification
# and save the data and plots under their own name
# We set the seed to easier be able to reproduce the result
UQ.quantify(figure_folder="figures_brunel_SR",
filename="brunel_SR",
seed=10)
We then change the parameters, and perform the uncertainty quantification and sensitivity analysis for the new set of parameters, again specifying a filename and figure folder.
# Change the set of parameters
UQ.parameters = parameters_AI
# Perform uncertainty quantification on the new parameter set
# and save the data and plots under their own name
# We set the seed to easier be able to reproduce the result
data = UQ.quantify(figure_folder="figures_brunel_AI",
filename="brunel_AI",
seed=10)
The complete code is:
import uncertainpy as un
import chaospy as cp
from brunel import brunel_network
# Create a Nest model from the brunel network function
# We set ``ignore=True`` since we are not interested in
# the model result itself.
# This is recommended for NEST models as long as you do not
# need the model results, since the uncertainty calculations for the
# for the model results require much time and memory.
model = un.NestModel(run=brunel_network, ignore=True)
# Parametes for the synchronous regular (SR) state
parameters = {"eta": cp.Uniform(1.5, 3.5),
"g": cp.Uniform(1, 3),
"delay": cp.Uniform(1.5, 3)}
parameters_SR = un.Parameters(parameters)
# Parameter for the asynchronous irregular (AI) state
parameters = {"eta": cp.Uniform(1.5, 2.2),
"g": cp.Uniform(5, 8),
"delay": cp.Uniform(1.5, 3)}
parameters_AI = un.Parameters(parameters)
# Initialize network features
features = un.NetworkFeatures()
# Set up the problem
UQ = un.UncertaintyQuantification(model,
parameters=parameters_SR,
features=features)
# Perform uncertainty quantification
# and save the data and plots under their own name
# We set the seed to easier be able to reproduce the result
UQ.quantify(figure_folder="figures_brunel_SR",
filename="brunel_SR",
seed=10)
# Change the set of parameters
UQ.parameters = parameters_AI
# Perform uncertainty quantification on the new parameter set
# and save the data and plots under their own name
# We set the seed to easier be able to reproduce the result
data = UQ.quantify(figure_folder="figures_brunel_AI",
filename="brunel_AI",
seed=10)