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

update multi-area-model.ipynb

parent 65fa0b19
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:b952d0ea tags:
#### Notebook structure <a class="anchor" id="toc"></a>
* [S0. Configuration](#section_0)
* [S1. Paramters specification](#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 connectivity](#section_2_3)
* [2.4. Run the simulation](#section_2_4)
* [S3. Data processing and simulation results analysis](#section_3)
* [S4. Simulation results visualization](#section_4)
%% Cell type:markdown id:bd3d4b0e tags:
<br>
%% Cell type:markdown id:d782e527 tags:
## S0. Configuration <a class="anchor" id="section_0"></a>
%% Cell type:code id:96517739 tags:
``` python
# Import dependencies
%matplotlib inline
import matplotlib.pyplot as plt
import numpy as np
import os
import nest
from IPython.display import display, HTML
# Import the MultiAreaModel class
from multiarea_model import MultiAreaModel
from config import base_path
```
%% Cell type:code id:2dd47c64 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:7e07b0d0 tags:
``` python
!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))
```
%% Cell type:markdown id:27160ba8 tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:565be233 tags:
<br>
%% Cell type:markdown id:df83f5ea-1c4b-44d3-9926-01786aa46e14 tags:
## S1. Paramters specification <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 |Variable |Value |Value description |
|:--------------------------------:|:--------------------------:|:------------------:|:------------------|
|Scaling factor |scale_down_to |0.005 |$^1$ |
|Ground state to metastable state |cc_weights_factor | |$^2$ |
|Replace non-simulated areas |replace_non_simulated_areas |'het_poisson_stat' |$^3$ |
|Simulation areas |sim_areas | |$^4$ |
%% Cell type:markdown id:a2161477 tags:
1. Scaling factor (`scale_down_to`) <br>
Scaling factor (scale_down_to) is the parameter 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.
Neurons and indegrees are both scaled down to 0.5%, 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. Ground state to metastable state (`cc_weights_factor`) <br>
Ground state to metastable state (cc_weights_factor) decides the switch of the model from ground state to metastable state. <br>
3. Replace non-simulated areas (`replace_non_simulated_areas`) <br>
Replace non-simulated areas (replace_non_simulated_areas) defines how non-simulated areas will be replaced. <br>
4. Simulation areas (`simulation_areas`) <br>
Simulation areas (simulation_areas) specifies the cortical areas included in the simulation process.
%% Cell type:code id:60265d52 tags:
``` python
# Model parameters
# Model paramters are most important among all the paramters, it directly affect the model itself and thus have a great impact on the simulation results. Model paramters define the connection, input, neuron, and network charateristics of the model, and therefore fall into four categories: **Connection paramters**, **Input paramters**, **Neuron paramters**, and **Network paramters**.
scale_down_to = 0.005 # Change it to 1. for running the fullscale network
cc_weights_factor =
# Scaling factor for cortico-cortical connections (chi)
# This scaling factor controls the cortico-cortical synaptic strength. If it is 1.0 then the inter-area synaptic strength is the same as the intra-areal. This factor changes the network activity from ground to metastable.
cc_weights_factor = 1.
replace_non_simulated_areas = 'het_poisson_stat'
sim_areas =
```
%% Cell type:markdown id:de11b07f tags:
### 1.2. Default parameters <a class="anchor" id="section_1_2"></a>
Default parameters ...<br>
We try out best not to confuse users with too many parameters. However, if you want to change the default parameters, 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)
conn_params = {'replace_non_simulated_areas': 'het_poisson_stat',
'g': -11.,
'K_stable': 'K_stable.npy',
'fac_nu_ext_TH': 1.2,
'fac_nu_ext_5E': 1.125,
'fac_nu_ext_6E': 1.41666667,
'av_indegree_V1': 3950.}
conn_params = {
# # It defines how non-simulated areas will be replaced
# Whether to replace non-simulated areas by Poisson sources with the same global rate rate_ext ('hom_poisson_stat') or by specific rates ('het_poisson_stat') or by time-varying specific current ('het_current_nonstat'). In the two latter cases, the data to replace the cortico-cortical input is loaded from `replace_cc_input_source`
'replace_non_simulated_areas': 'het_poisson_stat',
'g': -11.,
'K_stable': 'K_stable.npy',
# Increase the external input to 2/3E and 5E in area TH
'fac_nu_ext_TH': 1.2,
# Increase the external Poisson indegree onto 5E
'fac_nu_ext_5E': 1.125,
# Increase the external Poisson indegree onto 6E
'fac_nu_ext_6E': 1.41666667,
# Adjust the average indegree in V1 based on monkey data
'av_indegree_V1': 3950.
}
# Input parameters (input_params)
input_params = {'rate_ext': 10.}
# Neuron parameters (neuron_params)
neuron_params = {'V0_mean': -150.,
'V0_sd': 50.}
# Network parameters (network_params)
network_params = {'N_scaling': scale_down_to,
'K_scaling': scale_down_to,
'fullscale_rates': 'tests/fullscale_rates.json',
'input_params': input_params,
'connection_params': conn_params,
'neuron_params': neuron_params}
# Simulation parameters (sim_params)
sim_params = {'t_sim': 2000.,
'num_processes': 1,
'local_num_threads': 1,
'recording_dict': {'record_vm': False},
'rng_seed': 1} # global random seed
# Theory paramters (theory_params)
theory_params = {'dt': 0.1}
```
%% Cell type:markdown id:1472e9c5 tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:c532a861-824f-4713-a311-590aef8b6134 tags:
<br>
%% 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)
```
%% 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 "
"rate of {0:.3f} spikes/s across all populations.".format(np.mean(r[:, -1])))
```
%% Cell type:markdown id:2062ddf3 tags:
### 2.3. Extract 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>
**Warning**: memory explosion
%% Cell type:markdown id:b7396606 tags:
#### 2.3.1 Node indegrees
%% Cell type:code id:6316ac24 tags:
``` python
# Dictionary of nodes indegrees organized as:
# {<source_area>: {<source_pop>: {<target_area>: {<target_pop>: indegree_values}}}}
# M.K
```
%% Cell type:markdown id:253a2aba tags:
#### 2.3.2 Synapses
%% Cell type:code id:445a722a tags:
``` python
# 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:04894f5e-35ec-4b22-8891-bd7ba86098e9 tags:
<br>
%% Cell type:markdown id:0c1cad59-81d0-4e24-ac33-13c4ca8c6dec tags:
### 2.4. Run the 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()
```
%% Cell type:markdown id:fd6e3232 tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:4003c5a5-4a6f-49c5-be17-09f1bc68c411 tags:
<br>
%% Cell type:markdown id:28e071f8 tags:
## S3. Data processing and simulation results analysis <a class="anchor" id="section_3"></a>
%% Cell type:markdown id:89c7b7cf tags:
The following instructions will work when the `simulate` parameter is set to `True` during the creation of the MultiAreaModel object, and the `M.simulation.simulate()` method is executed.
%% Cell type:code id:dc3b1820 tags:
``` python
# Uncomment the lines in this code cell below to test if the number of synapses created by NEST matches the expected values
# """
# Test if the correct number of synapses has been created.
# """
# print("Testing synapse numbers")
# for target_area_name in M.area_list:
# target_area = M.simulation.areas[M.simulation.areas.index(target_area_name)]
# for source_area_name in M.area_list:
# source_area = M.simulation.areas[M.simulation.areas.index(source_area_name)]
# for target_pop in M.structure[target_area.name]:
# target_nodes = target_area.gids[target_pop]
# for source_pop in M.structure[source_area.name]:
# source_nodes = source_area.gids[source_pop]
# created_syn = nest.GetConnections(source=source_nodes,
# target=target_nodes)
# syn = M.synapses[target_area.name][target_pop][source_area.name][source_pop]
# assert(len(created_syn) == int(syn))
```
%% Cell type:markdown id:57401110 tags:
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:8726a93d tags:
### 1. Load spike data
%% Cell type:code id:cb8e3edd tags:
``` python
data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + "-spikes-1-0.dat", skiprows=3)
```
%% Cell type:markdown id:8793e033 tags:
### 2. Compute instantaneous rate per neuron across all populations
%% Cell type:code id:9590223b tags:
``` python
tsteps, spikecount = np.unique(data[:,1], return_counts=True)
rate = spikecount / M.simulation.params['dt'] * 1e3 / np.sum(M.N_vec)
```
%% Cell type:markdown id:ca603daf tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:b1320ab1 tags:
Go back to [Notebook structure](#toc)
%% Cell type:markdown id:529b1ade tags:
<br>
%% Cell type:markdown id:57ff902c-d6ce-4f96-9e4f-8e3e7166ab66 tags:
## S4. Simulation results visualization <a class="anchor" id="section_4"></a>
%% Cell type:markdown id:38ddd973 tags:
### 4.1. Instantaneous and mean rate
%% Cell type:code id:bea30fc8 tags:
``` python
fig, ax = plt.subplots()
ax.plot(tsteps, rate)
ax.plot(tsteps, np.average(rate)*np.ones(len(tsteps)), label='mean')
ax.set_title('instantaneous rate across all populations')
ax.set_xlabel('time (ms)')
ax.set_ylabel('rate (spikes / s)')
ax.set_xlim(0, sim_params['t_sim'])
ax.set_ylim(0, 50)
ax.legend()
```
%% 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