Skip to content
Snippets Groups Projects
Commit efad086c authored by Didi Hou's avatar Didi Hou
Browse files

update multi-area-model.ipynb

parent 346cccd6
No related branches found
No related tags found
1 merge request!35Pre-release MAM v1.1.0
%% Cell type:markdown id:b1331599 tags:
# Down-scaled multi-area model
%% Cell type:markdown id:edec8345-aec1-419e-b9e3-7f612aff8262 tags:
<img src="model_construction.png" alt="Model overview" width="1000"/>
%% Cell type:markdown id:f4a649cc-3b68-49e4-b2b6-6f29f13a6d9c tags:
The code in this notebook implements the down-scaled version of spiking network model of macaque visual cortex developed at the Institute of Neuroscience and Medicine (INM-6), Research Center Jülich. The full-scale model has been documented in the following publications:
1. Schmidt M, Bakker R, Hilgetag CC, Diesmann M & van Albada SJ
Multi-scale account of the network structure of macaque visual cortex
Brain Structure and Function (2018), 223: 1409 [https://doi.org/10.1007/s00429-017-1554-4](https://doi.org/10.1007/s00429-017-1554-4)
2. Schuecker J, Schmidt M, van Albada SJ, Diesmann M & Helias M (2017)
Fundamental Activity Constraints Lead to Specific Interpretations of the Connectome.
PLOS Computational Biology, 13(2): e1005179. [https://doi.org/10.1371/journal.pcbi.1005179](https://doi.org/10.1371/journal.pcbi.1005179)
3. Schmidt M, Bakker R, Shen K, Bezgin B, Diesmann M & van Albada SJ (2018)
A multi-scale layer-resolved spiking network model of
resting-state dynamics in macaque cortex. PLOS Computational Biology, 14(9): e1006359. [https://doi.org/10.1371/journal.pcbi.1006359](https://doi.org/10.1371/journal.pcbi.1006359)
%% Cell type:markdown id:ff3671a8 tags:
<br>
%% 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. Instantiate a multi-area model](#section_2_1)
* [2.2. Predict firing rates from theory](#section_2_2)
* [2.3. Extract and visualize interareal connectivity](#section_2_3)
* [2.4. Run a simulation](#section_2_4)
* [S3. Simulation Results Visualization](#section_3)
* [3.1. Instantaneous and mean firing rate across all populations](#section_3_1)
* [3.2. Resting state plots](#section_3_2)
* [3.3. Time-averaged population rates](#section_3_3)
%% Cell type:markdown id:620501a2 tags:
<br>
%% 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 numpy as np
import os
import nest
import json
import sys
from IPython.display import display, HTML
import warnings
sys.path.append('./figures/MAM2EBRAINS')
from multiarea_model import MultiAreaModel
from multiarea_model import Analysis
from M2E_visualize_interareal_connectivity import visualize_interareal_connectivity
from M2E_visualize_instantaneous_and_mean_firing_rates import plot_instan_mean_firing_rate
from M2E_visualize_resting_state import plot_resting_state
from config import base_path, data_path
```
%% Output
-- N E S T --
Copyright (C) 2004 The NEST Initiative
Version: 3.5
Built: Jul 12 2023 06:25:27
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
%%capture captured
!pip install nested_dict dicthash
```
%% Cell type:code id:1d440c07-9b69-4e52-8573-26b13493bc5a tags:
``` python
# Jupyter notebook display format setting
style = """
<style>
table {float:left}
</style>
"""
display(HTML(style))
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] |0.005 |$^1$ |
|cc_weights_factor|1. |[1.0, 2.5] |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
# Value range/options: (0, 1.]
# Value assgined: 0.005
scale_down_to = 0.005 # Change it to 1. for running the fullscale network
# Scaling factor for cortico-cortical connections (chi)
# Value range/options: [1., 2.5]
# Value assgined: 1.0
cc_weights_factor = 1.0
# Cortical areas included in the simulation
# Value range/options: any sublist of complete_ares_list
# Value assgined: complete_ares_list
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
# Value range/options: None, 'hom_poisson_stat', 'het_poisson_stat', 'het_current_nonstat'
# Value assgined: 'het_poisson_stat'
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
# 't_sim': 1500., # 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. Instantiate a multi-area model <a class="anchor" id="section_2_1"></a>
%% Cell type:code id:ab25f9f8 tags:
``` python
# %%capture captured
M = MultiAreaModel(network_params,
simulation=True,
sim_spec=sim_params,
theory=True,
theory_spec=theory_params)
```
%% Output
Initializing network from dictionary.
RAND_DATA_LABEL 5313
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])))
```
%% Cell type:markdown id:2062ddf3 tags:
### 2.3. Extract and visualize 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
# Neuron numbers
# Dictionary of neuron numbers
# print(M.N)
# Array of neuron numbers
# (M.N_vec)
# print(M.N_vec)
```
%% Cell type:code id:8408d463-557b-481b-afc1-5fbbbd67306d tags:
``` python
# Indegrees
# Dictionary of nodes indegrees organized as:
# {<source_area>: {<source_pop>: {<target_area>: {<target_pop>: indegree_values}}}}
# M.K
# print(M.K)
# Array of nodes indegrees
# M.K_matrix.shape
# print(M.K_matrix.shape)
```
%% 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
# print(M.synapses)
# Array of
# M.syn_matrix
# print(M.syn_matrix)
```
%% Cell type:code id:05512922-26e5-425f-90a4-0df7c2279ccf tags:
``` python
visualize_interareal_connectivity(M)
```
%% Output
Initializing network from dictionary.
RAND_DATA_LABEL 3491
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
--------------------
{}
========================================
Initializing network from dictionary.
RAND_DATA_LABEL 5677
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
--------------------
{}
========================================
%% Cell type:markdown id:bae85d86-157c-47a2-9826-860b410a440e tags:
Comparable figure in our publications: <br>
1. Schmidt M, Bakker R, Hilgetag CC, Diesmann M & van Albada SJ
Multi-scale account of the network structure of macaque visual cortex
Brain Structure and Function (2018), 223: 1409 [https://doi.org/10.1007/s00429-017-1554-4](https://doi.org/10.1007/s00429-017-1554-4)
**Fig. 4 D Area-level connectivity of the model, based on data in a–c, expressed as relative indegrees for each target area**
%% 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
# %%capture captured
# 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.
Sep 14 14:49:07 SimulationManager::set_status [Info]:
Temporal resolution changed from 0.1 to 0.1 ms.
Rank 0: created area V1 with 0 local nodes
Memory after V1 : 1519.98 MB
Rank 0: created area V2 with 0 local nodes
Memory after V2 : 1546.70 MB
Rank 0: created area VP with 0 local nodes
Memory after VP : 1575.62 MB
Rank 0: created area V3 with 0 local nodes
Memory after V3 : 1603.86 MB
Rank 0: created area V3A with 0 local nodes
Memory after V3A : 1623.84 MB
Rank 0: created area MT with 0 local nodes
Memory after MT : 1649.37 MB
Rank 0: created area V4t with 0 local nodes
Memory after V4t : 1672.66 MB
Rank 0: created area V4 with 0 local nodes
Memory after V4 : 1699.76 MB
Rank 0: created area VOT with 0 local nodes
Memory after VOT : 1725.05 MB
Rank 0: created area MSTd with 0 local nodes
Memory after MSTd : 1746.52 MB
Rank 0: created area PIP with 0 local nodes
Memory after PIP : 1767.86 MB
Rank 0: created area PO with 0 local nodes
Memory after PO : 1789.35 MB
Rank 0: created area DP with 0 local nodes
Memory after DP : 1809.61 MB
Rank 0: created area MIP with 0 local nodes
Memory after MIP : 1831.09 MB
Rank 0: created area MDP with 0 local nodes
Memory after MDP : 1852.60 MB
Rank 0: created area VIP with 0 local nodes
Memory after VIP : 1874.41 MB
Rank 0: created area LIP with 0 local nodes
Memory after LIP : 1898.46 MB
Rank 0: created area PITv with 0 local nodes
Memory after PITv : 1923.76 MB
Rank 0: created area PITd with 0 local nodes
Memory after PITd : 1948.95 MB
Rank 0: created area MSTl with 0 local nodes
Memory after MSTl : 1970.45 MB
Rank 0: created area CITv with 0 local nodes
Memory after CITv : 1989.50 MB
Rank 0: created area CITd with 0 local nodes
Memory after CITd : 2008.86 MB
Rank 0: created area FEF with 0 local nodes
Memory after FEF : 2030.36 MB
Rank 0: created area TF with 0 local nodes
Memory after TF : 2045.99 MB
Rank 0: created area AITv with 0 local nodes
Memory after AITv : 2068.57 MB
Rank 0: created area FST with 0 local nodes
Memory after FST : 2085.26 MB
Rank 0: created area 7a with 0 local nodes
Memory after 7a : 2106.58 MB
Rank 0: created area STPp with 0 local nodes
Memory after STPp : 2125.29 MB
Rank 0: created area STPa with 0 local nodes
Memory after STPa : 2144.34 MB
Rank 0: created area 46 with 0 local nodes
Memory after 46 : 2159.67 MB
Rank 0: created area AITd with 0 local nodes
Memory after AITd : 2182.33 MB
Rank 0: created area TH with 0 local nodes
Memory after TH : 2195.03 MB
Created areas and internal connections in 2.22 seconds.
Created cortico-cortical connections in 22.01 seconds.
Sep 14 14:49:31 NodeManager::prepare_nodes [Info]:
Preparing 20780 nodes for simulation.
Sep 14 14:49:31 SimulationManager::start_updating_ [Info]:
Number of local nodes: 20780
Simulation time (ms): 2000
Number of OpenMP threads: 1
Number of MPI processes: 1
Sep 14 14:51:01 SimulationManager::run [Info]:
Simulation finished.
Simulated network in 89.83 seconds.
%% Cell type:markdown id:fd6e3232 tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:bb71c922 tags:
## S3. Simulation Results Visualziation <a class="anchor" id="section_3"></a>
%% Cell type:code id:c1d7aa61-e85a-4e6a-9e01-e018413a572b tags:
``` python
# Instantiate an analysis class and load spike data
A = Analysis(network=M,
simulation=M.simulation,
data_list=['spikes'],
load_areas=None)
```
%% Output
loading spikes
%% Cell type:markdown id:38ddd973 tags:
### 3.1. Instantaneous and mean firing rate across all populations <a class="anchor" id="section_3_1"></a>
%% Cell type:code id:bea30fc8 tags:
``` python
plot_instan_mean_firing_rate(M)
```
%% Output
%% Cell type:markdown id:e91c436e-db94-4cd7-a531-29c032efeeae tags:
### 3.2 Resting state plots <a class="anchor" id="section_3_2"></a>
%% Cell type:markdown id:aeae56a4 tags:
Comparable figure in our publications: <br>
3. Schmidt M, Bakker R, Shen K, Bezgin B, Diesmann M & van Albada SJ (2018)
A multi-scale layer-resolved spiking network model of
resting-state dynamics in macaque cortex. PLOS Computational Biology, 14(9): e1006359. [https://doi.org/10.1371/journal.pcbi.1006359](https://doi.org/10.1371/journal.pcbi.1006359)
**Fig 5. Resting state of the model with χ =1.9.**
%% Cell type:code id:ae19bcc3 tags:
``` python
# Choose 3 areas from the complete_area_list to show their sipking activity
# By default, it's set as ['V1', 'V2', 'FEF']
raster_areas = ['V2', 'V1', 'FEF']
plot_resting_state(M, A, data_path, raster_areas)
```
%% Output
Loading data from file
Loading data from file
Computing rate time series done
Loading data from file
Computing synchrony done
pop_LvR
pop_rates
synchrony
%% Cell type:markdown id:473d0882-8e45-4330-bfa2-2c7e1af0dac4 tags:
### 3.3 Time-averaged population rates <a class="anchor" id="section_4_3"></a>
An 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
# %%capture captured
A.show_rates()
```
%% Output
%% Cell type:markdown id:b03d44e8-2216-44ff-ada4-83e9c3e6d30a tags:
|Area index|0|1|2|3|4|5|6|7|8|9|10|11|12|13|14|15|16|17|18|19|20|21|22|23|24|25|26|27|28|29|30|31|
|:--------:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|Area |V1|V2|VP|V3|PIP|V3A|MT|V4t|V4|PO|VOT|DP|MIP|MDP|MSTd|VIP|LIP|PITv|PITd|AITv|MSTl|FST|CITv|CITd|7a|STPp|STPa|FEF|46|TF|TH|AITd|
%% Cell type:markdown id:ef74ca3e-98dc-49c9-a4a0-2c640e29b1d9 tags:
Go back to [Notebook structure](#toc)
......
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