Skip to content
Snippets Groups Projects
Commit 9862c79d authored by Didi Hou's avatar Didi Hou Committed by Administrator
Browse files

/

parent 2e6ebb2b
No related branches found
No related tags found
1 merge request!35Pre-release MAM v1.1.0
import numpy as np
def load_data(M):
# load spike data and calculate instantaneous and mean firing rates
data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + "-spikes-1-0.dat", skiprows=3)
......
import numpy as np
def load_data(M):
# load spike data and calculate instantaneous and mean firing rates
data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + "-spikes-1-0.dat", skiprows=3)
......
%% Cell type:markdown id:b1331599 tags:
# Down-scaled multi-area model
%% Cell type:markdown id:b952d0ea tags:
#### Notebook structure <a class="anchor" id="toc"></a>
* [S0. Configuration](#section_0)
* [S1. Parameterization](#section_1)
* [1.1. Parameters to tune](#section_1_1)
* [1.2. Default parameters](#section_1_2)
* [S2. Multi-Area Model Instantiation and Simulation](#section_2)
* [2.1. Insantiate a multi-area model](#section_2_1)
* [2.2. Predict firing rates from theory](#section_2_2)
* [2.3. Extract interareal connectivity](#section_2_3)
* [2.4. Run a simulation](#section_2_4)
* [S3. SExtract Interneural Connectivity](#section_3)
* [S4. Data Loading and Processing](#section_4)
* [S5. Simulation Results Visualziation](#section_5)
* [5.1. Instantaneous and mean firing rate across all populations](#section_5_1)
* [5.2 Resting state plots](#section_5_2)
* [5.3 Time-averaged population rates](#section_5_3)
%% Cell type:markdown id:d782e527 tags:
## S0. Configuration <a class="anchor" id="section_0"></a>
%% Cell type:code id:9d6cc7d9-3110-4d96-9f9a-9ec7dee6d145 tags:
``` python
# Create config file
with open('config.py', 'w') as fp:
fp.write(
'''import os
base_path = os.path.abspath(".")
data_path = os.path.abspath("simulations")
jobscript_template = "python {base_path}/run_simulation.py {label}"
submit_cmd = "bash -c"
''')
```
%% Cell type:code id:96517739 tags:
``` python
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import os
import nest
import json
from multiarea_model import MultiAreaModel
from multiarea_model import Analysis
from config import base_path, data_path
import sys
sys.path.append('./figures')
```
%% Output
-- N E S T --
Copyright (C) 2004 The NEST Initiative
Version: 3.4
Built: May 17 2023 20:48:31
This program is provided AS IS and comes with
NO WARRANTY. See the file LICENSE for details.
Problems or suggestions?
Visit https://www.nest-simulator.org
Type 'nest.help()' to find out more about NEST.
%% Cell type:code id:7e07b0d0 tags:
``` python
!pip install nested_dict dicthash
```
%% Output
Requirement already satisfied: nested_dict in /srv/main-spack-instance-2302/spack/var/spack/environments/ebrains-23-02/.spack-env/._view/6axslmv6jvf4v2nte3uwlayg4vhsjoha/lib/python3.8/site-packages (1.61)
Requirement already satisfied: dicthash in /srv/main-spack-instance-2302/spack/var/spack/environments/ebrains-23-02/.spack-env/._view/6axslmv6jvf4v2nte3uwlayg4vhsjoha/lib/python3.8/site-packages (0.0.2)
Requirement already satisfied: future in /srv/main-spack-instance-2302/spack/var/spack/environments/ebrains-23-02/.spack-env/._view/6axslmv6jvf4v2nte3uwlayg4vhsjoha/lib/python3.8/site-packages (from dicthash) (0.18.2)
%% Cell type:code id:1d440c07-9b69-4e52-8573-26b13493bc5a tags:
``` python
# Jupyter notebook display format setting
from IPython.display import display, HTML
style = """
<style>
table {float:left}
</style>
"""
display(HTML(style))
# Ignore and don't display warning messages
import warnings
warnings.filterwarnings("ignore")
```
%% Output
%% Cell type:markdown id:27160ba8 tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:df83f5ea-1c4b-44d3-9926-01786aa46e14 tags:
## S1. Parameterization <a class="anchor" id="section_1"></a>
%% Cell type:markdown id:30655817 tags:
### 1.1. Parameters to tune <a class="anchor" id="section_1_1"></a>
%% Cell type:markdown id:4f67c1ba tags:
|Parameter |Default value |Value range/options |Value assigned |Description |
|:----------------------------:|:-----------------------:|:--------------------------------------------------------------------:|:------------------:|:-----------:|
|scale_down_to |1. |(0, 1.] |0.005 |$^1$ |
|cc_weights_factor |1. |(0, 1.] |1. |$^2$ |
|areas_simulated |complete_area_list |Sublists of complete_area_list |complete_area_list |$^3$ |
|replace_non_simulated_areas |None |None, 'hom_poisson_stat', 'het_poisson_stat', 'het_current_nonstat' |'het_poisson_stat' |$^4$ |
%% Cell type:markdown id:a2161477 tags:
1. `scale_down_to` <br>
`scale_down_to` is the down-scaling factor which defines the the ratio of the full scale multi-area model being down-scaled to a model with fewer neurons and indegrees so as to be simulated on machines with lower computational ability and the simulation results can be obtained within relative shorter period of time. <br> Its deafualt value if `1.` meaning full scale simulation. <br> In the pre-set downscale version, it's set as `0.005`, where the numer of neurons and indegrees are both scaled down to 0.5% of its full scale amount, where the model can usually be simulated on a local machine. <br> **Warning**: This will not yield reasonable dynamical results from the network and is only meant to demonstrate the simulation workflow <br>
2. `cc_weights_factor` <br>
This scaling factor controls the cortico-cortical synaptic strength. <br> By default it's set as `1.0`, where the inter-area synaptic strength is the same as the intra-areal. <br> **Important**: This factor changes the network activity from ground state to metastable state. <br>
3. `areas_simulated` <br>
This parameter specifies the cortical areas included in the simulation process. Its default value is `complete_area_list` meaning all the areas in the complete_area_list will be actually simulated. <br>
complete_area_list = ['V1', 'V2', 'VP', 'V3', 'V3A', 'MT', 'V4t', 'V4', 'VOT', 'MSTd', 'PIP', 'PO', 'DP', 'MIP', 'MDP', 'VIP', 'LIP', 'PITv', 'PITd', 'MSTl', 'CITv', 'CITd', 'FEF', 'TF', 'AITv', 'FST', '7a', 'STPp', 'STPa', '46', 'AITd', 'TH'] <br>
The value assigned to simulation_areas can be any sublist of the compete_area_list specifying areas a user want to include in his/her simulation. <br>
4. `replace_non_simulated_areas` <br>
The paramter `replace_non_simulated_areas` defines how non-simulated areas will be replaced. <br> It's set as `None` by default when the parameter areas_simulated is set as full_area_list where all areas will be simulated so that no areas need to be replaced. <br> Other options are: `'hom_poisson_stat'`, `'het_poisson_stat'`, and `'het_current_nonstat'`. `'hom_poisson_stat'` is a manually set parameter which can be tuned. When it's set as 'het_poisson_stat' or 'het_current_nonstat', the data to replace the cortico-cortical input is loaded from 'replace_cc_input_source' which is the firing rates of our full scale simulation results. The differenc between 'het_poisson_stat' and 'het_current_nonstat' is that 'het_poisson_stat' is the mean of the time-series firing rate so that it's static, yet 'het_current_nonstat' is time-varying specific current, which is varying by time.
%% Cell type:code id:60265d52 tags:
``` python
# Downscaling factor
scale_down_to = 0.005 # Change it to 1. for running the fullscale network
# Scaling factor for cortico-cortical connections (chi)
cc_weights_factor = 1.
# Cortical areas included in the simulation
areas_simulated = ['V1', 'V2', 'VP', 'V3', 'V3A', 'MT', 'V4t', 'V4', 'VOT', 'MSTd', 'PIP', 'PO', 'DP', 'MIP', 'MDP', 'VIP', 'LIP', 'PITv', 'PITd', 'MSTl', 'CITv', 'CITd', 'FEF', 'TF', 'AITv', 'FST', '7a', 'STPp', 'STPa', '46', 'AITd', 'TH']
# Firing rates used to replace the non-simulated areas
replace_non_simulated_areas = 'het_poisson_stat'
```
%% Cell type:markdown id:de11b07f tags:
### 1.2. Default parameters <a class="anchor" id="section_1_2"></a>
We try our best not to confuse users with too many parameters. However, if you want to change more parameters and explore the model, you can do so by passing a dictionary to the `default_params` argument of the `MultiAreaModel` class.
%% Cell type:code id:6e4bed8d tags:
``` python
# Connection parameters
conn_params = {
'replace_non_simulated_areas': 'het_poisson_stat', # Whether to replace non-simulated areas by Poisson sources with the same global rate, by default: None
'g': -11., # It sets the relative inhibitory synaptic strength, by default: -16.
'K_stable': 'K_stable.npy', # Whether to apply the stabilization method of Schuecker, Schmidt et al. (2017), by default: None
'fac_nu_ext_TH': 1.2, # Increase the external input to 2/3E and 5E in area TH
'fac_nu_ext_5E': 1.125, # Increase the external Poisson indegree onto 5E
'fac_nu_ext_6E': 1.41666667, # Increase the external Poisson indegree onto 6E
'av_indegree_V1': 3950. # Adjust the average indegree in V1 based on monkey data
}
# Input parameters
input_params = {
'rate_ext': 10. # Rate of the Poissonian spike generator (in spikes/s)
}
# Neuron parameters
neuron_params = {
'V0_mean': -150., # Mean for the distribution of initial membrane potentials, by default: -100.
'V0_sd': 50.} # Standard deviation for the distribution of initial membrane potentials, by default: 50.
# Network parameters
network_params = {
'N_scaling': scale_down_to, # Scaling of population sizes, by default: 1.
'K_scaling': scale_down_to, # Scaling of indegrees, by default: 1.
'fullscale_rates': 'tests/fullscale_rates.json', # Absolute path to the file holding full-scale rates for scaling synaptic weights, by default: None
'input_params': input_params, # Input parameters
'connection_params': conn_params, # Connection parameters
'neuron_params': neuron_params # Neuron parameters
}
# Simulation parameters
sim_params = {
'areas_simulated': areas_simulated,
't_sim': 2000., # Simulated time (in ms), by default: 10.0
'num_processes': 1, # The number of MPI processes, by default: 1
'local_num_threads': 1, # The number of threads per MPI process, by default: 1
'recording_dict': {'record_vm': False},
'rng_seed': 1 # global random seed
}
# Theory paramters (theory_params)
theory_params = {
'dt': 0.1 # The time step of the mean-field theory integration, by default: 0.01
}
```
%% Cell type:markdown id:1472e9c5 tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:de4a6703 tags:
## S2. Multi-Area Model Instantiation and Simulation <a class="anchor" id="section_2"></a>
%% Cell type:markdown id:1fd58841 tags:
### 2.1. Insantiate a multi-area model <a class="anchor" id="section_2_1"></a>
%% Cell type:code id:ab25f9f8 tags:
``` python
M = MultiAreaModel(network_params,
simulation=True,
sim_spec=sim_params,
theory=True,
theory_spec=theory_params)
```
%% Output
Initializing network from dictionary.
RAND_DATA_LABEL 1109
RAND_DATA_LABEL 8335
Error in library("aod") : there is no package called ‘aod’
Execution halted
No R installation or IndexError, taking hard-coded SLN fit parameters.
========================================
Customized parameters
--------------------
{'K_scaling': 0.005,
'N_scaling': 0.005,
'connection_params': {'K_stable': 'K_stable.npy',
'av_indegree_V1': 3950.0,
'fac_nu_ext_5E': 1.125,
'fac_nu_ext_6E': 1.41666667,
'fac_nu_ext_TH': 1.2,
'g': -11.0,
'replace_non_simulated_areas': 'het_poisson_stat'},
'fullscale_rates': 'tests/fullscale_rates.json',
'input_params': {'rate_ext': 10.0},
'neuron_params': {'V0_mean': -150.0, 'V0_sd': 50.0}}
========================================
Simulation label: 27d81076e6d6e9e591684be053078477
Copied files.
Initialized simulation class.
%% Cell type:markdown id:91649c30 tags:
### 2.2. Predict firing rates from theory <a class="anchor" id="section_2_2"></a>
%% Cell type:code id:6a7ddf0e tags:
``` python
p, r = M.theory.integrate_siegert()
print("Mean-field theory predicts an average "
"firing rate of {0:.3f} spikes/s across all populations.".format(np.mean(r[:, -1])))
```
%% Output
Iteration: 0
Mean-field theory predicts an average firing rate of 29.588 spikes/s across all populations.
%% Cell type:markdown id:2062ddf3 tags:
### 2.3. Extract interareal connectivity <a class="anchor" id="section_2_3"></a>
%% Cell type:markdown id:8a7c09e0 tags:
The connectivity and neuron numbers are stored in the attributes of the model class. Neuron numbers are stored in `M.N` as a dictionary (and in `M.N_vec` as an array), indegrees in `M.K` as a dictionary (and in `M.K_matrix` as an array). Number of synapses can also be access via `M.synapses` (and in `M.syn_matrix` as an array). <br>
%% Cell type:code id:6316ac24 tags:
``` python
# Indegrees
# Dictionary of nodes indegrees organized as:
# {<source_area>: {<source_pop>: {<target_area>: {<target_pop>: indegree_values}}}}
# M.K
```
%% Cell type:code id:445a722a tags:
``` python
# Synapses
# Dictionary of synapses that target neurons receive, it is organized as:
# {<source_area>: {<source_pop>: {<target_area>: {<target_pop>: number_of_synapses}}}}
# M.synapses
```
%% Cell type:markdown id:e67f37e9-ec8d-4bb1-bd21-45e966f47ab6 tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:0c1cad59-81d0-4e24-ac33-13c4ca8c6dec tags:
### 2.4. Run a simulation <a class="anchor" id="section_2_4"></a>
%% Cell type:code id:15778e9c tags:
``` python
# run the simulation, depending on the model parameter and downscale ratio, the running time varies largely.
M.simulation.simulate()
```
%% Output
Prepared simulation in 0.00 seconds.
Rank 0: created area V1 with 0 local nodes
Memory after V1 : 1911.96 MB
Rank 0: created area V2 with 0 local nodes
Memory after V2 : 1938.56 MB
Rank 0: created area VP with 0 local nodes
Memory after VP : 1967.84 MB
Memory after VP : 1967.79 MB
Rank 0: created area V3 with 0 local nodes
Memory after V3 : 1996.08 MB
Rank 0: created area V3A with 0 local nodes
Memory after V3A : 2016.02 MB
Rank 0: created area MT with 0 local nodes
Memory after MT : 2041.53 MB
Rank 0: created area V4t with 0 local nodes
Memory after V4t : 2066.47 MB
Memory after V4t : 2066.51 MB
Rank 0: created area V4 with 0 local nodes
Memory after V4 : 2093.50 MB
Memory after V4 : 2093.58 MB
Rank 0: created area VOT with 0 local nodes
Memory after VOT : 2118.69 MB
Memory after VOT : 2118.77 MB
Rank 0: created area MSTd with 0 local nodes
Memory after MSTd : 2140.32 MB
Memory after MSTd : 2140.24 MB
Rank 0: created area PIP with 0 local nodes
Memory after PIP : 2161.67 MB
Memory after PIP : 2161.59 MB
Rank 0: created area PO with 0 local nodes
Memory after PO : 2183.14 MB
Rank 0: created area DP with 0 local nodes
Memory after DP : 2203.41 MB
Rank 0: created area MIP with 0 local nodes
Memory after MIP : 2224.87 MB
Memory after MIP : 2224.91 MB
Rank 0: created area MDP with 0 local nodes
Memory after MDP : 2246.34 MB
Memory after MDP : 2246.42 MB
Rank 0: created area VIP with 0 local nodes
Memory after VIP : 2268.28 MB
Memory after VIP : 2268.23 MB
Rank 0: created area LIP with 0 local nodes
Memory after LIP : 2292.38 MB
Memory after LIP : 2292.30 MB
Rank 0: created area PITv with 0 local nodes
Memory after PITv : 2317.55 MB
Memory after PITv : 2317.61 MB
Rank 0: created area PITd with 0 local nodes
Memory after PITd : 2342.85 MB
Memory after PITd : 2342.81 MB
Rank 0: created area MSTl with 0 local nodes
Memory after MSTl : 2364.20 MB
Memory after MSTl : 2364.27 MB
Rank 0: created area CITv with 0 local nodes
Memory after CITv : 2383.38 MB
Memory after CITv : 2383.45 MB
Rank 0: created area CITd with 0 local nodes
Memory after CITd : 2402.79 MB
Rank 0: created area FEF with 0 local nodes
Memory after FEF : 2424.18 MB
Memory after FEF : 2424.29 MB
Rank 0: created area TF with 0 local nodes
Memory after TF : 2439.86 MB
Memory after TF : 2439.82 MB
Rank 0: created area AITv with 0 local nodes
Memory after AITv : 2462.58 MB
Memory after AITv : 2462.50 MB
Rank 0: created area FST with 0 local nodes
Memory after FST : 2479.19 MB
Memory after FST : 2479.23 MB
Rank 0: created area 7a with 0 local nodes
Memory after 7a : 2500.48 MB
Memory after 7a : 2500.40 MB
Rank 0: created area STPp with 0 local nodes
Memory after STPp : 2519.23 MB
Rank 0: created area STPa with 0 local nodes
Memory after STPa : 2538.38 MB
Memory after STPa : 2538.30 MB
Rank 0: created area 46 with 0 local nodes
Memory after 46 : 2553.75 MB
Rank 0: created area AITd with 0 local nodes
Memory after AITd : 2576.30 MB
Memory after AITd : 2576.34 MB
Rank 0: created area TH with 0 local nodes
Memory after TH : 2588.98 MB
Created areas and internal connections in 2.29 seconds.
Created cortico-cortical connections in 22.97 seconds.
Simulated network in 73.53 seconds.
Memory after TH : 2589.02 MB
Created areas and internal connections in 2.21 seconds.
Created cortico-cortical connections in 22.15 seconds.
Simulated network in 78.61 seconds.
%% Cell type:markdown id:fd6e3232 tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:28e071f8 tags:
## S3. Extract Interneural Connectivity <a class="anchor" id="section_3"></a>
%% Cell type:markdown id:57401110 tags:
**Warning**: Memory explosion <br>
To obtain the connections information, you can extract the lists of connected sources and targets. Moreover, you can access additional synaptic details, such as synaptic weights and delays.
%% Cell type:code id:e7eb052e tags:
``` python
# conns = nest.GetConnections()
# conns_sparse_matrix = conns.get(['source', 'target', 'weight'])
# srcs = conns_sparse_matrix['source']
# tgts = conns_sparse_matrix['target']
# weights = conns_sparse_matrix['weight']
```
%% Cell type:markdown id:ef4b2e4b tags:
You can determine the area and subpopulation to which the neuron ID ranges belong by referring to the file `network_gids.txt`, which is automatically generated during network creation.
%% Cell type:code id:902f2800 tags:
``` python
# # Open the file using a with statement
# with open(os.path.join(M.simulation.data_dir,"recordings/network_gids.txt"), "r") as file:
# # Read the contents of the file
# gids = file.read()
# # Print the contents
# print(gids)
```
%% Cell type:markdown id:b1320ab1 tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:57ff902c-d6ce-4f96-9e4f-8e3e7166ab66 tags:
## S4. Data Loading and Processing <a class="anchor" id="section_4"></a>
%% Cell type:code id:f5b58845-4d1a-430f-83f4-402fdf918aef tags:
``` python
label_spikes = M.simulation.label
label = M.simulation.label
from MAM2EBRAINS_LOAD_DATA import load_data
A, tsteps, firing_rate = load_data(M)
```
%% Output
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
NameError Traceback (most recent call last)
Input In [14], in <module>
2 label = M.simulation.label
4 from MAM2EBRAINS_LOAD_DATA import load_data
----> 5 A, tsteps, firing_rate = load_data(M)
TypeError: load_data() missing 1 required positional argument: 'A'
File ~/MAM2EBRAINS/./figures/MAM2EBRAINS_LOAD_DATA.py:3, in load_data(M)
1 def load_data(M):
2 # load spike data and calculate instantaneous and mean firing rates
----> 3 data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + "-spikes-1-0.dat", skiprows=3)
4 tsteps, spikecount = np.unique(data[:,1], return_counts=True)
5 firing_rate = spikecount / M.simulation.params['dt'] * 1e3 / np.sum(M.N_vec)
NameError: name 'np' is not defined
%% Cell type:markdown id:2da9728d-4481-4a15-b810-d125e39cbe4e tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:bb71c922 tags:
## S5. Simulation Results Visualziation <a class="anchor" id="section_5"></a>
%% Cell type:markdown id:38ddd973 tags:
### 5.1. Instantaneous and mean firing rate across all populations <a class="anchor" id="section_5_1"></a>
%% Cell type:code id:bea30fc8 tags:
``` python
from MAM2EBRAINS_VISUALIZATION import plot_instan_mean_firing_rate
plot_instan_mean_firing_rate(tsteps, firing_rate, sim_params)
```
%% Cell type:markdown id:e91c436e-db94-4cd7-a531-29c032efeeae tags:
### 5.2 Resting state plots <a class="anchor" id="section_5_2"></a>
%% Cell type:markdown id:aeae56a4 tags:
**Fig 5. Resting state of the model with χ =1.9.** (A-C) Raster plot of spiking activity of 3% of the neurons in area V1 (A), V2 (B), and FEF (C). Blue: excitatory neurons, red: inhibitory neurons. (D-F) Spiking statistics across all 32 areas for the respective populations shown as area-averaged box plots. Crosses: medians, boxes: interquartile range (IQR), whiskers extend to the most extremeobservat ions within 1.5×IQR beyond the IQR. (D) Population-averaged firing rates. (E) Average pairwise correlation coefficients of spiking activity. (F) Irregularity measured by revised local variation LvR averaged across neurons. (G) Area-averaged firing rates, shown as raw binned spike histograms with 1ms bin width (gray) and convolved histograms, with aGaussian kernel (black) of optimal width.
%% Cell type:code id:ae19bcc3 tags:
``` python
from MAM2EBRAINS_VISUALIZATION import plot_resting_state
plot_resting_state(A, label_spikes)
```
%% Cell type:markdown id:473d0882-8e45-4330-bfa2-2c7e1af0dac4 tags:
### 5.3 Time-averaged population rates <a class="anchor" id="section_5_3"></a>
Plot overview over time-averaged population rates encoded in colors with areas along x-axis and populations along y-axis.
%% Cell type:code id:721d1f03-df25-468d-8075-a807025a9c58 tags:
``` python
# area_list = ['V1', 'V2', 'VP', 'V3', 'V3A', 'MT', 'V4t', 'V4', 'VOT', 'MSTd', 'PIP', 'PO', 'DP', 'MIP', 'MDP', 'VIP', 'LIP', 'PITv', 'PITd', 'MSTl', 'CITv', 'CITd', 'FEF', 'TF', 'AITv', 'FST', '7a', 'STPp', 'STPa', '46', 'AITd', 'TH']
# output = {'pdf', 'png', 'eps'}, optional
A.show_rates(area_list)
```
%% Cell type:markdown id:ef74ca3e-98dc-49c9-a4a0-2c640e29b1d9 tags:
Go back to [Notebook structure](#toc)
......
source diff could not be displayed: it is too large. Options to address this: view the blob.
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment