multi-area-model.ipynb

 
 Down-scaled multi-area model
 
 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:

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
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
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
 Notebook Outline 

S0. Configuration
S1. Parameterization
1.1. Parameters to tune
S2. Model Configuration, Instantiation and Simulation
2.1. Configuring model parameters
2.2. Instantiate a multi-area model
2.3. Predict firing rates from theory
2.4. Extract and visualize inter-areal connectivity
2.5. Run a simulation
S3. Visualization of Network Dynamics 
3.1. Mean firing rate over simulated populations
3.2. Instantaneous firing rate over simulated areas
3.3. Time-averaged firing rate over all populations
3.4. Network dynamics
3.5. Functional connectivity
Additional Notes
 S0. Configuration 
 # 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"
''')
 
 %matplotlib inline
import numpy as np
import sys
import os
from IPython.display import display, HTML
import warnings

from multiarea_model import MultiAreaModel
from multiarea_model import Analysis
from config import base_path, data_path

sys.path.append('./figures/MAM2EBRAINS')
 
 # Jupyter notebook display format setting
style = """
<style>
table {float:left}
</style>
"""
display(HTML(style))

warnings.filterwarnings('ignore')
 
 Go back to Notebook Outline
 S1. Parameterization 
 1.1. Parameters to tune 
 The values assigned for the following parameters are kept the same as in the paper except for the scale_down_to which is set as 0.006 enabling to simulate a down-scaled multi-area model with 2GB RAM. By default, it is set to 1.0 for simulating the full-scale model.
 

Parameter
Default value
Value range/options
Value assigned
Description


scale_down_to
$1.0$
$(0, 1.0]$
$0.006$
$^1$


cc_weights_factor
$1.9$
$\geq 1.0$
$1.9$
$^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$


g
$-11.0$
$\leq -1.0$
$-11.0$
$^5$


rate_ext
$10.0$
$\geq 0.0$
$10.0$
$^6$


scale_down_to is the down-scaling factor that defines the ratio by which the full-scale multi-area model is reduced to a model with fewer neurons and indegrees. This reduction is essential to enable simulation on machines with limited computational power, ensuring that simulation results can be obtained in a relatively shorter timeframe. 
 If the value is scale_down_to = 1., the full-scale network will be simulated. 
 In the pre-set down-scaled version, scale_down_to = 0.006. This setting reduces the number of neurons and indegrees to 0.6 % of their full-scale counterparts, facilitating simulation on a typical local machine. 
 Warning: This may not yield reasonable results from the network dynamics and is only meant to demonstrate the simulation workflow! 
 
cc_weights_factor is the scaling factor that controls the cortico-cortical synaptic strength. 
 By default it is set to 1.9, keeping the same value for producing the metastable state as in the original paper. 
 Important: This factor plays a crucial role in transitioning the network activity from the ground to the metastable state.  In the full-scale network, the ground state and metastable state activities are achieved when this parameter is set to 1.0 and 1.9, respectively. In the down-scaled multi-area model, a similar metastable state may not be achieved or achieved with a different value. 
areas_simulated specifies the cortical areas to be included in the simulation process. Its default value is complete_area_list meaning all the areas in the complete_area_list will be simulated.
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']

The value assigned to areas_simulated can be any sublist of complete_area_list.

replace_non_simulated_areas defines how non-simulated areas will be replaced. 
 When all areas are included, it is set as None by default. 
 Other options are: 'hom_poisson_stat', 'het_poisson_stat', and 'het_current_nonstat'.
 'hom_poisson_stat' replaces the non-simulated areas by Poisson sources with the same global rate rate_ext. The 'het_poisson_stat' and 'het_current_nonstat' options use the loaded specific rates from 'replace_cc_input_source', which contains the area-specific firing rates of our full-scale simulation results. The difference is that 'het_poisson_stat' replaces the non-simulated areas by Poisson spike trains and 'het_current_nonstat' replaces it with a time-varying current input.
g defines the relative inhibitory synaptic strength (in relative units to the excitatory synaptic strength). By default: -11.0, as used in the full-scale network. g = -1.0 means equal excitatory and inhibitory strengths, and g < -1.0 results in stronger inhibition than excitation.
rate_ext defines the rate of the Poissonian spike generator (in spikes/s), by default: 10.0. It also serves as one of the input parameters of the model. When a larger value is assigned to rate_ext, the excitatory background noise is increased. Note that the external Poisson indegree onto 5E and 6E is increased by a factor of 1.125 and 1.41666667 respectively, and the external Poisson indegree onto 23E and 5E in area TH is increased by a factor of 1.2.
 # Downscaling factor
# value range/options: (0, 1.], change it to 1. to simulate the full-scale network
scale_down_to = 0.006

# Scaling factor for cortico-cortical connections (Chi) 
# value range/options: [1., 2.5], 
# a weight factor of 1.0 produces Ground state activity.
# 1.9 was assigned to produce results in Schmidt et al. (2018).
cc_weights_factor = 1.9

# Cortical areas included in the simulation
# value range/options: any sublist of complete_area_list
# where complete_area_list is
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']
areas_simulated = complete_area_list
# areas_simulated = ['V1', 'V2']

# Firing rates used to replace the non-simulated areas
# value range/options: None, 'hom_poisson_stat', 'het_poisson_stat', 'het_current_nonstat'
# if areas_simulated is complete_area_list, then replace_non_simulated_areas will be set as None 
# regardless of the value assigned below
replace_non_simulated_areas = 'het_poisson_stat'

# Relative inhibitory synaptic strength (in relative units), by default: -11.
g = -11.

# Rate of the Poissonian spike generator (in spikes/s), by default: 10.
rate_ext = 10.
 
 Go back to Notebook Outline
 S2. Model Configuration, Instantiation and Simulation 
 2.1. Configuring model parameters 
 We try our best not to confuse users with too many parameters. So, the few parameters tunned will be automatically assigned in this section to properly run the simulation. 
However, if you want to explore the model, you can alter other parameters related to the network or simulation configuration by passing them in the network_params and sim_params dictionaries below. If this is not the case, you can execute the cell the way it is.
 # Determine replace_cc_input_source
replace_cc_input_source = None                                               # By default, replace_cc_input_source is set to None
                                                                             # where areas_simulated is complete_area_list                                                           
if set(areas_simulated) != set(complete_area_list):                                                                                       
    if replace_non_simulated_areas == 'hom_poisson_stat':                   
        replace_cc_input_source = None
    elif replace_non_simulated_areas == 'het_poisson_stat' or replace_non_simulated_areas == 'het_current_nonstat':
        replace_cc_input_source = os.path.join(base_path, 'tests/fullscale_rates.json')
    else:
        raise Exception("'hom_poisson_stat', 'het_poisson_stat', or 'het_current_nonstat' should be assigned to replace_non_simulated_areas when not all areas are simulated!")

# Determine cc_weights_I_factor from cc_weights_factor
if cc_weights_factor == 1.0:                                                  # For ground state with cc_weights_factor as 1., 
    cc_weights_I_factor = 1.0                                                 # cc_weights_I_factor is set to 1.
elif cc_weights_factor > 1.0:                                                 # For cc_weights_factor larger than 1.0,
    cc_weights_I_factor = 2.0                                                 # cc_weights_I_factor is set to 2.
else:                                                                         # cc_weights_factor outside of (1., 2.5], raise error
    raise Exception("A value that is equal to or larger than 1.0 should be assigned to the parameter cc_weights_factor!")

# Connection parameters
conn_params = {
    'replace_non_simulated_areas': replace_non_simulated_areas,               # Whether to replace non-simulated areas by Poisson sources 
    'g': g,                                                                   # It sets the relative inhibitory synaptic strength, by default: -11.
    'replace_cc_input_source': replace_cc_input_source,                       # Specify the data used to replace non-simulated areas      
    'cc_weights_factor': cc_weights_factor,
    'cc_weights_I_factor': cc_weights_I_factor
}

# Input parameters
input_params = {
    'rate_ext': rate_ext                                                      # Rate of the Poissonian spike generator (in spikes/s), by default: 10.
} 

# Network parameters
network_params = {
    'N_scaling': scale_down_to,                                               # Scaling of population sizes, by default: 1. for full scale multi-area model
    'K_scaling': scale_down_to,                                               # Scaling of indegrees, by default: 1. for full scale multi-area model
    'fullscale_rates': os.path.join(base_path, '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
} 

# Simulation parameters
sim_params = {
    'areas_simulated': areas_simulated,                                       # Cortical areas included in the simulation
    't_sim': 2000.,                                                           # Simulated time (in ms), by default: 10.
    'rng_seed': 1                                                             # Global random seed
}
 
 Go back to Notebook Outline
 2.2. Instantiate a multi-area model 
 M = MultiAreaModel(network_params, 
                   simulation=True,
                   sim_spec=sim_params,
                   theory=True)
 
 2.3. Predict firing rates from theory 
 Note: the prediction may differ from the simulation results, especially in the presence of synchrony.
 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])))
 
 Go back to Notebook Outline
 2.4. Extract and visualize inter-areal connectivity 
 The connectivity and neuron numbers are stored in the attributes of the model class. 

Neuron numbers of all populations in each area are stored in M.N as a dictionary (and in M.N_vec as an array).
Indegrees are stored in M.K as a dictionary (and in M.K_matrix as an array).

Dictionary of nodes indegrees organized as:

{<target_area>: {<target_pop>: {<source_area>: {<source_pop>: indegree_values}}}}
Number of synapses can be accessed via M.synapses (and in M.syn_matrix as an array). 

Dictionary of synapses that target neurons receive, it is organized as:

{<target_area>: {<target_pop>: {<source_area>: {<source_pop>: number_of_synapses}}}}
 The figure below shows the inter-areal connectivity of the model expressed as the relative indegrees of each target area. The relative indegree of a target area from a specific source area is calculated by dividing its indegree by the sum of indegrees that the target area receives from all sources.
 Comparable figure in our publications: 

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 

Fig. 4D Area-level connectivity of the model, based on data in a–c, expressed as relative indegrees for each target area.
 from M2E_visualize_interareal_connectivity import visualize_interareal_connectivity
visualize_interareal_connectivity(M)
 
 Go back to Notebook Outline
 2.5. Run a simulation 
 # Run the simulation, depending on the model parameter and downscale ratio, the running time varies largely.
M.simulation.simulate()
 
 Reminder: The spike trains of simulated results are saved to the folder with path ./simulations/<simulation_label>/recordings where the <simulation_label> is displayed in the output of 2.2. All statistics describing network dynamics are computed from the saved spike trains.
 Go back to Notebook Outline
 S3. Visualization of Network Dynamics 
 Important: cc_weights_factor plays a crucial role in transitioning the network activity from the ground to the metastable state.  In the full-scale network, the ground state and metastable state activities are achieved when this parameter is set to 1.0 and 1.9, respectively.
 # Create an instance of Analysis to load data
A = Analysis(M, M.simulation, data_list=['spikes'], load_areas=None)
 
 3.1. Mean firing rate over simulated populations 
 # Print the mean firing rate over simulated populations
from M2E_firing_rate import mean_firing_rate
mean_firing_rate(M, data_path)
 
 3.2. Instantaneous firing rate over simulated areas 
 from M2E_firing_rate import plot_firing_rate_over_areas
plot_firing_rate_over_areas(M, data_path)
 
 3.3. Time-averaged firing rate over simulated populations 
 An overview of time-averaged firing rate over simulated populations encoded in colors with areas along x-axis and populations along y-axis. The cells of population 4E and 4I in area TH are labeled with X as area TH does not have layer 4.
 from M2E_visualize_time_ave_pop_rates import plot_time_averaged_population_rates
plot_time_averaged_population_rates(M, data_path)
 
 3.4. Network dynamics 
 Comparable figures in our publications: 

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 

Fig 3.  Ground state of the model. 

Fig 5.  Resting state of the model with χ = 1.9 (metastable state).
 # Choose at most 3 areas from the areas_simulated to show their spiking activities
# By default, the list is ['V1', 'V2', 'FEF'] when all areas from complete_area_list are simulated
raster_areas = ['V1', 'V2', 'FEF']

from M2E_visualize_dynamics import visual_dynamics
visual_dynamics(M, data_path, raster_areas)
 
 3.5. Functional connectivity 
 Comparison of area-level functional connectivity (FC) between the down-scaled MAM and macaque experimental data. (A) Simulated FC measured by the zero-time-lag correlation coefficient of synaptic input currents. (B) FC of macaque resting-state fMRI (see Materials and methods).
 from M2E_visualize_fc import visualize_fc
visualize_fc(M, data_path)
 
 Go back to Notebook Outline
 Additional Notes 
 
Simulation data 

The spike data of all simulated populations for all simulations are saved in ./simulations/<simulation_label>/recordings where <simulation_label> can be accessed in the output of 2.2. Or users can see their latest simulation by checking the column "Last Modified" and find the folder with the latest change.
Statistics 

The statistics of network dynamics computed from the spike trains can be found in ./simulations/<simulation_label>/Analysis. You may also find more statistics defined in ./multiarea_model/analysis.py to explore more about network dynamics.
Scripts for visualizing network dynamics 

The scripts for computing statistics and plotting the figures in S3 can be found in ./figures/MAM2EBRAINS.
 Go back to Notebook Outline