Skip to content
Snippets Groups Projects
Commit cbf48d7c authored by didihou's avatar didihou Committed by Administrator
Browse files

/

parent 45151de9
No related branches found
No related tags found
1 merge request!35Pre-release MAM v1.1.0
Expand all Hide whitespace changes
Inline Side-by-side
.ipynb_checkpoints/multi-area-model-checkpoint.ipynb
+ 16
− 16
View file @ cbf48d7c
This diff is collapsed.
figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_resting_state-checkpoint.py
+ 18
− 13
View file @ cbf48d7c
......@@ -17,6 +17,8 @@ from matplotlib import gridspec
icolor = myred
ecolor = myblue
from M2E_LOAD_DATA import load_and_create_data
def set_boxplot_props(d):
for i in range(len(d['boxes'])):
if i % 2 == 0:
......@@ -32,8 +34,11 @@ def set_boxplot_props(d):
pl.setp(d['means'], marker='x', color='k',
markerfacecolor='k', markeredgecolor='k', markersize=3.)
def plot_resting_state(M, A, data_path):
# label_spikes = M.simulation.label
def plot_resting_state(M, data_path):
# load data
# A = load_and_create_data(M)
label_spikes = M.simulation.label
label = M.simulation.label
t_sim = M.simulation.params["t_sim"]
......@@ -148,17 +153,17 @@ def plot_resting_state(M, A, data_path):
# """
# M = MultiAreaModel({})
# spike data
# spike_data = {}
# for area in areas:
# spike_data[area] = {}
# for pop in M.structure[area]:
# spike_data[area][pop] = np.load(os.path.join(data_path,
# label_spikes,
# 'recordings',
# '{}-spikes-{}-{}.npy'.format(label_spikes,
# area, pop)))
spike_data = A.spike_data
spike data
spike_data = {}
for area in areas:
spike_data[area] = {}
for pop in M.structure[area]:
spike_data[area][pop] = np.load(os.path.join(data_path,
label_spikes,
'recordings',
'{}-spikes-{}-{}.npy'.format(label_spikes,
area, pop)))
# spike_data = A.spike_data
# stationary firing rates
fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json')
......
figures/MAM2EBRAINS/M2E_visualize_resting_state.py
+ 18
− 13
View file @ cbf48d7c
......@@ -17,6 +17,8 @@ from matplotlib import gridspec
icolor = myred
ecolor = myblue
from M2E_LOAD_DATA import load_and_create_data
def set_boxplot_props(d):
for i in range(len(d['boxes'])):
if i % 2 == 0:
......@@ -32,8 +34,11 @@ def set_boxplot_props(d):
pl.setp(d['means'], marker='x', color='k',
markerfacecolor='k', markeredgecolor='k', markersize=3.)
def plot_resting_state(M, A, data_path):
# label_spikes = M.simulation.label
def plot_resting_state(M, data_path):
# load data
# A = load_and_create_data(M)
label_spikes = M.simulation.label
label = M.simulation.label
t_sim = M.simulation.params["t_sim"]
......@@ -148,17 +153,17 @@ def plot_resting_state(M, A, data_path):
# """
# M = MultiAreaModel({})
# spike data
# spike_data = {}
# for area in areas:
# spike_data[area] = {}
# for pop in M.structure[area]:
# spike_data[area][pop] = np.load(os.path.join(data_path,
# label_spikes,
# 'recordings',
# '{}-spikes-{}-{}.npy'.format(label_spikes,
# area, pop)))
spike_data = A.spike_data
spike data
spike_data = {}
for area in areas:
spike_data[area] = {}
for pop in M.structure[area]:
spike_data[area][pop] = np.load(os.path.join(data_path,
label_spikes,
'recordings',
'{}-spikes-{}-{}.npy'.format(label_spikes,
area, pop)))
# spike_data = A.spike_data
# stationary firing rates
fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json')
......
multi-area-model.ipynb
+ 5
− 6
View file @ cbf48d7c
%% 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 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: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
from multiarea_model import MultiAreaModel
from config import base_path, data_path
sys.path.append('./figures/MAM2EBRAINS')
```
%% 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))
```
%% 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. Insantiate 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
Error in library("aod") : there is no package called ‘aod’
Execution halted
%% 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 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
# M.N
# Array of neuron numbers
# 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
# Array of nodes indegrees
# 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
# Array of
# M.syn_matrix
```
%% Cell type:code id:05512922-26e5-425f-90a4-0df7c2279ccf tags:
``` python
%%capture captured
# from M2E_visualize_interareal_connectivity import visualize_interareal_connectivity
# visualize_interareal_connectivity(M)
```
%% Cell type:markdown id:bae85d86-157c-47a2-9826-860b410a440e tags:
Based on: <br>
Schmidt M, Bakker R, Hilgetag CC, Diesmann M & van Albada SJ <br>
Multi-scale account of the network structure of macaque visual cortex <br>
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) <br>
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()
```
%% 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: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
from M2E_visualize_instantaneous_and_mean_firing_rates import plot_instan_mean_firing_rate
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:
**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 M2E_visualize_resting_state import plot_resting_state
plot_resting_state(M, A, label_spikes, data_path)
plot_resting_state(M, data_path)
```
%% Output
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [15], line 1
----> 1 from M2E_visualize_resting_state import plot_resting_state
2 plot_resting_state(M, A, label_spikes, data_path)
File ~/MAM2EBRAINS/./figures/MAM2EBRAINS/M2E_visualize_resting_state.py:21
18 ecolor = myblue
20 # label_spikes = M.simulation.label
---> 21 label = M.simulation.label
23 def set_boxplot_props(d):
24 for i in range(len(d['boxes'])):
NameError: name 'M' is not defined
Cell In [16], line 2
1 from M2E_visualize_resting_state import plot_resting_state
----> 2 plot_resting_state(M, A, label_spikes, data_path)
NameError: name 'A' is not defined
%% 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()
```
%% Cell type:code id:5b40db5d-51b2-4d16-b36b-9f1995452b05 tags:
``` python
|Index|0|1|2|3|4|5|6|7|
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|
```
%% Cell type:markdown id:b03d44e8-2216-44ff-ada4-83e9c3e6d30a tags:
|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