Merge pull request #62 from didi-hou/plot_enhance

Plots enhancing

Merge pull request #62 from didi-hou/plot_enhance
Plots enhancing
af205934 · shimoura · GitHub · 681447b3 · cb9c0fe5 · af205934
Unverified Commit af205934 authored 1 year ago by shimoura Committed by GitHub 1 year ago
--- a/figures/MAM2EBRAINS/M2E_firing_rate.py
+++ b/figures/MAM2EBRAINS/M2E_firing_rate.py
@@ -70,18 +70,18 @@ def plot_firing_rate_over_areas(M, data_path):
        matrix.append(binned_spikes)
    matrix = np.array(matrix)
-    normalized_matrix = (matrix - np.min(matrix)) / (np.max(matrix) - np.min(matrix))
+    # normalized_matrix = (matrix - np.min(matrix)) / (np.max(matrix) - np.min(matrix))
-    fig = pl.figure()
+    fig = pl.figure(figsize=(12, 5))
-    fig.suptitle('Normalized instantaneous firing rate over simulated areas', fontsize=16, x=0.45, y=0.95)
+    fig.suptitle('Instantaneous firing rate over simulated areas', fontsize=16, x=0.45, y=0.95)
    ax = pl.subplot()
    cmap = pl.get_cmap('YlOrBr')
    # masked_matrix = np.ma.masked_where(np.isnan(matrix), matrix)
    ax.patch.set_hatch('x')
-    im = ax.pcolormesh(normalized_matrix, cmap=cmap, edgecolors='None')
+    im = ax.pcolormesh(matrix, cmap=cmap, edgecolors='None', vmin=0)
-    ax.set_xlim(0, normalized_matrix[0].size)
+    ax.set_xlim(0, matrix[0].size)
    ax.set_xticks([i for i in np.arange(0, 100, 10)])
    ax.set_xticklabels([i for i in np.arange(500, 600, 10)])

--- a/figures/MAM2EBRAINS/M2E_visualize_interareal_connectivity.py
+++ b/figures/MAM2EBRAINS/M2E_visualize_interareal_connectivity.py
@@ -28,7 +28,7 @@ def visualize_interareal_connectivity(M):
    nrows = 1
    ncols = 2
    # width = 6.8556
-    width = 15
+    width = 12
    panel_wh_ratio = 0.7 * (1. + np.sqrt(5)) / 2.  # golden ratio
    height = width / panel_wh_ratio * float(nrows) / ncols
@@ -36,13 +36,13 @@ def visualize_interareal_connectivity(M):
    pl.rcParams['figure.figsize'] = (width, height)
    fig = pl.figure()
-    fig.suptitle('Interareal probabilistic connectivity for full-scale (left) and down-scaled (right) multi-area model', fontsize=16, y=1.05)
+    fig.suptitle('Area-level connectivity of the full-scale and down-scaled MAM expressed as relative indegrees for each target area', fontsize=15, x=0.5, y=1.05)
    axes = {}
    # gs1 = gridspec.GridSpec(2, 2)
    gs1 = gridspec.GridSpec(1, 2)
    # gs1.update(left=0.06, right=0.95, top=0.95, bottom=0.1, wspace=0.1, hspace=0.3)
-    gs1.update(left=0.06, right=0.95, top=0.95, bottom=0.1, wspace=0.3, hspace=0.3)
+    gs1.update(left=0.1, right=0.95, top=0.95, bottom=0.1, wspace=0.3, hspace=0.3)
    # axes['A'] = pl.subplot(gs1[:1, :1])
    # axes['B'] = pl.subplot(gs1[:1, 1:2])

--- a/figures/MAM2EBRAINS/M2E_visualize_time_ave_pop_rates.py
+++ b/figures/MAM2EBRAINS/M2E_visualize_time_ave_pop_rates.py
@@ -29,7 +29,7 @@ def plot_time_averaged_population_rates(M, data_path, area_list=None, **keywords
    # matrix = np.zeros((len(area_list), len(A.network.structure['V1'])))
    matrix = np.zeros((len(area_list), len(M.structure['V1'])))
-    fig = plt.figure(figsize=(12, 4))
+    fig = plt.figure(figsize=(12, 3))
    fig.suptitle('Time-averaged firing rate over simulated populations', fontsize=15, x=0.43)
    ax = fig.add_subplot(111)
@@ -51,20 +51,23 @@ def plot_time_averaged_population_rates(M, data_path, area_list=None, **keywords
                rate = np.nan
            matrix[i][j] = rate
-    cm = plt.cm.jet
+    # cm = plt.cm.jet
-    cm = cm.from_list('mycmap', [(0., 64./255., 192./255.),  # custom dark blue
+    # cm = cm.from_list('mycmap', [(0., 64./255., 192./255.),  # custom dark blue
-                                 (0., 128./255., 192./255.),  # custom light blue
+    #                              (0., 128./255., 192./255.),  # custom light blue
-                                 'white',
+    #                              'white',
-                                 (245./255., 157./255., 115./255.),  # custom light red
+    #                              (245./255., 157./255., 115./255.),  # custom light red
-                                 (192./255., 64./255., 0.)], N=256)  # custom dark red
+    #                              (192./255., 64./255., 0.)], N=256)  # custom dark red
+    cm = plt.get_cmap('YlOrBr')
    cm.set_under('0.3')
-    cm.set_bad('k')
+    # cm.set_bad('k')
+    cm.set_bad('white')
    matrix = np.transpose(matrix)
    masked_matrix = np.ma.masked_where(np.isnan(matrix), matrix)
    ax.patch.set_hatch('x')
-    im = ax.pcolormesh(masked_matrix, cmap=cm, edgecolors='None', norm=LogNorm(
+    # im = ax.pcolormesh(masked_matrix, cmap=cm, edgecolors='None', norm=LogNorm(
-        vmin=0.01, vmax=100.))
+    #     vmin=0.01, vmax=100.))
+    im = ax.pcolormesh(masked_matrix, cmap=cm, edgecolors='None', vmin=0)
    ax.set_xlim(0, matrix[0].size)
    x_index = np.arange(4.5, 31.6, 5.0)
@@ -82,9 +85,17 @@ def plot_time_averaged_population_rates(M, data_path, area_list=None, **keywords
    ax.set_yticklabels(['6I', '6E', '5I', '5E', '4I', '4E', '2/3I', '2/3E'])
    ax.set_ylabel('Population', size=13)
    ax.set_xlabel('Area index', size=13)
-    t = FixedLocator([0.01, 0.1, 1., 10., 100.])
+    # t = FixedLocator([0.01, 0.1, 1., 10., 100.])
+    # t = FixedLocator([0, 10, 20, 30, 40, 50])
+    # Iterate over the data and add 'X' for masked cells
+    for i in range(masked_matrix.shape[0]):
+        for j in range(masked_matrix.shape[1]):
+            if masked_matrix.mask[i, j]:
+                ax.text(j + 0.5, i + 0.5, 'X', va='center', ha='center', color='black', fontsize=23)
-    plt.colorbar(im, ticks=t)
+    # plt.colorbar(im, ticks=t)
+    plt.colorbar(im)
    if 'output' in keywords:
        # plt.savefig(os.path.join(A.output_dir, '{}_rates.{}'.format(A.simulation.label,

--- a/multi-area-model.ipynb
+++ b/multi-area-model.ipynb
@@ -537,7 +537,9 @@
  {
   "cell_type": "markdown",
   "id": "2714fd2b-df6e-45b9-a7f4-f240a8c65ecf",
-   "metadata": {},
+   "metadata": {
+    "tags": []
+   },
   "source": [
    "### 3.2. Normalized instantaneous firing rate over simulated areas <a class=\"anchor\" id=\"section_3_2\"></a>"
   ]
@@ -566,7 +568,7 @@
   "id": "b6746ea4-91bd-44f5-acb6-83df15b05480",
   "metadata": {},
   "source": [
-    "An overview of time-averaged firing rate over simulated populations encoded in colors with areas along x-axis and populations along y-axis."
+    "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."
   ]
  },
  {

 %% 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:b952d0ea tags:
 ## Notebook Outline <a class="anchor" id="toc"></a>
 * [S0. Configuration](#section_0)
 * [S1. Parameterization](#section_1)
    * [1.1. Parameters to tune](#section_1_1)
 * [S2. Model Configuration, Instantiation and Simulation](#section_2)
    * [2.1. Configuring model parameters](#section_2_1)
    * [2.2. Instantiate a multi-area model](#section_2_2)
    * [2.3. Predict firing rates from theory](#section_2_3)
    * [2.4. Extract and visualize interareal connectivity](#section_2_3)
    * [2.5. Run a simulation](#section_2_5)
 * [S3. Visualization of Network Dynamics](#section_3)
    * [3.1. Mean firing rate over simulated populations](#section_3_1)
    * [3.2. Normalized instantaneous firing rate over simulated areas](#section_3_2)
    * [3.3. Time-averaged firing rate over all populations](#section_3_3)
    * [3.4. Network dynamics](#section_3_4)
    * [3.5. Functional connectivity](#section_3_5)
 * [Additional Notes](#section_4)
 %% 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
 !pip install python-louvain
 !pip install community
 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')
 ```
 %% Cell type:code id:06764b48-a3b0-4518-ba92-816398bb72b0 tags:
 ``` python
 # Jupyter notebook display format setting
 style = """
 <style>
 table {float:left}
 </style>
 """
 display(HTML(style))
 warnings.filterwarnings('ignore')
 ```
 %% Cell type:markdown id:27160ba8 tags:
 Go back to [Notebook Outline](#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$        |$(0, 1.0]$         |$0.006$       |$^1$       |
 |cc_weights_factor|$1.0$        |$[1.0, 2.5]$       |$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$ |
 %% Cell type:markdown id:a2161477 tags:
 1. `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. <br> If the value is `scale_down_to = 1.`, the full-scale network will be simulated. <br> 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. <br> **Warning**: This may not yield reasonable results from the network dynamics and is only meant to demonstrate the simulation workflow! <br>
 2. `cc_weights_factor` is the scaling factor that controls the cortico-cortical synaptic strength. <br> By default it is set to `1.0`, ensuring that inter- and intra-area synapses have equal strength. <br> **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.<br>
 3. `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.
 ```python
 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.
 4. `replace_non_simulated_areas` defines how non-simulated areas will be replaced. <br> When all areas are included, it is set as `None` by default. <br> Other options are: `'hom_poisson_stat'`, `'het_poisson_stat'`, and `'het_current_nonstat'`.<br> `'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.
 5. `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.
 6. `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.
 %% Cell type:code id:60265d52 tags:
 ``` python
 # Downscaling factor
 # Value range/options: (0, 1.], change it to 1. to run 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.
 # A value larger than 1. reproduces the "Metastable" 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.
 ```
 %% Cell type:markdown id:1472e9c5 tags:
 Go back to [Notebook Outline](#toc)
 %% Cell type:markdown id:de4a6703 tags:
 ## S2. Model Configuration, Instantiation and Simulation <a class="anchor" id="section_2"></a>
 %% Cell type:markdown id:9608b6d9-c7e4-4b2d-9c2b-0d43c6415b48 tags:
 ### 2.1. Configuring model parameters <a class="anchor" id="section_2_1"></a>
 %% Cell type:markdown id:56e2b6f6-d67a-4a74-973b-c8507a975bf8 tags:
 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.
 %% Cell type:code id:fc49ba27-cb8e-441a-bf34-432971555b90 tags:
 ``` python
 # 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 and cc_weights_factor <= 2.5:                    # For metastable state with cc_weights_factor in (1., 2.5],
    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 in range [1.0, 2.5] 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
 }
 ```
 %% Cell type:markdown id:b3b5f634-284d-4b7b-88f6-17d26fb7743c tags:
 Go back to [Notebook Outline](#toc)
 %% Cell type:markdown id:1fd58841 tags:
 ### 2.2. Instantiate a multi-area model <a class="anchor" id="section_2_2"></a>
 %% Cell type:code id:ab25f9f8 tags:
 ``` python
 M = MultiAreaModel(network_params,
                   simulation=True,
                   sim_spec=sim_params,
                   theory=True)
 ```
 %% Cell type:markdown id:91649c30 tags:
 ### 2.3. Predict firing rates from theory <a class="anchor" id="section_2_3"></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:c361fa92-4f50-4519-9592-60a39888a12b tags:
 Go back to [Notebook Outline](#toc)
 %% Cell type:markdown id:2062ddf3 tags:
 ### 2.4. Extract and visualize interareal connectivity <a class="anchor" id="section_2_4"></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 are stored in `M.K` as a dictionary (and in `M.K_matrix` as an array).<br>
  Dictionary of nodes indegrees organized as:<br>
  `{<source_area>: {<source_pop>: {<target_area>: {<target_pop>: indegree_values}}}}`
 - Number of synapses can be accessed via `M.synapses` (and in `M.syn_matrix` as an array). <br>
  Dictionary of synapses that target neurons receive, it is organized as:<br>
  `{<source_area>: {<source_pop>: {<target_area>: {<target_pop>: number_of_synapses}}}}`
 %% Cell type:markdown id:5c15ee30-b5eb-4024-b127-98ab68337ec0 tags:
 Comparable figure in our publications: <br>
 1. Schmidt M, Bakker R, Hilgetag CC, Diesmann M & van Albada SJ <br>
   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) <br>
   **Fig. 4D Area-level connectivity of the model, based on data in a–c, expressed as relative indegrees for each target area.**
 %% Cell type:code id:05512922-26e5-425f-90a4-0df7c2279ccf tags:
 ``` python
 from M2E_visualize_interareal_connectivity import visualize_interareal_connectivity
 visualize_interareal_connectivity(M)
 ```
 %% Cell type:markdown id:e67f37e9-ec8d-4bb1-bd21-45e966f47ab6 tags:
 Go back to [Notebook Outline](#toc)
 %% Cell type:markdown id:0c1cad59-81d0-4e24-ac33-13c4ca8c6dec tags:
 ### 2.5. Run a simulation <a class="anchor" id="section_2_5"></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:72a89aeb-78b9-4642-be0d-6de4d6c8f817 tags:
 **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.
 %% Cell type:markdown id:fd6e3232 tags:
 Go back to [Notebook Outline](#toc)
 %% Cell type:markdown id:bb71c922 tags:
 ## S3. Visualization of Network Dynamics <a class="anchor" id="section_3"></a>
 %% Cell type:markdown id:5d150e0d-9174-4737-be1d-0fa2ac282419 tags:
 **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.
 %% Cell type:code id:44473002-b970-4f48-ac3d-1e01093133fd tags:
 ``` python
 # Create an instance of Analysis to load data
 A = Analysis(M, M.simulation, data_list=['spikes'], load_areas=None)
 ```
 %% Cell type:markdown id:38ddd973 tags:
 ### 3.1. Mean firing rate over simulated populations <a class="anchor" id="section_3_1"></a>
 %% Cell type:code id:bea30fc8 tags:
 ``` python
 # Print the mean firng rate over simulated populations
 from M2E_firing_rate import mean_firing_rate
 mean_firing_rate(M, data_path)
 ```
 %% Cell type:markdown id:2714fd2b-df6e-45b9-a7f4-f240a8c65ecf tags:
 ### 3.2. Normalized instantaneous firing rate over simulated areas <a class="anchor" id="section_3_2"></a>
 %% Cell type:code id:61eb9d0d-06ef-454e-98c3-645152208979 tags:
 ``` python
 from M2E_firing_rate import plot_firing_rate_over_areas
 plot_firing_rate_over_areas(M, data_path)
 ```
 %% Cell type:markdown id:de0317c8-1fcd-405c-b0c9-36336aa1f3ad tags:
 ### 3.3. Time-averaged firing rate over simulated populations <a class="anchor" id="section_3_3"></a>
 %% Cell type:markdown id:b6746ea4-91bd-44f5-acb6-83df15b05480 tags:
-An overview of time-averaged firing rate over simulated populations encoded in colors with areas along x-axis and populations along y-axis.
+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.
 %% Cell type:code id:14fa641f-0371-413e-945c-c8d4579aead6 tags:
 ``` python
 from M2E_visualize_time_ave_pop_rates import plot_time_averaged_population_rates
 plot_time_averaged_population_rates(M, data_path)
 ```
 %% Cell type:markdown id:e91c436e-db94-4cd7-a531-29c032efeeae tags:
 ### 3.4. Network dynamics <a class="anchor" id="section_3_4"></a>
 %% Cell type:markdown id:2ab1e9e5-3eb8-40b2-a4e9-2e749978878d tags:
 Comparable figures in our publications: <br>
 1. 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) <br>
   **Fig 3.  Ground state of the model.** <br>
   **Fig 5.  Resting state of the model with χ = 1.9 (metastable state).**
 %% Cell type:code id:ae19bcc3 tags:
 ``` python
 # Choose at most 3 areas from the areas_simulated to show their sipking activities
 # By default, it's set as ['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)
 ```
 %% Cell type:markdown id:b53058b5-c0bd-4837-8289-9226872317cc tags:
 ### 3.5. Functional connectivity <a class="anchor" id="section_3_5"></a>
 %% Cell type:markdown id:013adaf8-af8b-470e-94f0-b69121d1ca2c tags:
 Comparison of inter-area interactions of the down-scaled MAM visualized as functional connectivity and experimental data. (A) Simulated functional connectivity (FC) for w = 1.9 measured by the zero-time-lag correlation coefficient of synaptic input currents. (B) FC of macaque resting-state fMRI (see Materials and methods).
 %% Cell type:markdown id:678741b2-6dfa-4f1d-b4d8-ccd84f104767 tags:
 **Order of cortical areas**: <br>
 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
 %% Cell type:code id:b4b25621-618d-4594-8cf7-9b9002837d69 tags:
 ``` python
 from M2E_visualize_fc import visualize_fc
 label = M.simulation.label
 visualize_fc(M, data_path, label)
 ```
 %% Cell type:markdown id:ef74ca3e-98dc-49c9-a4a0-2c640e29b1d9 tags:
 Go back to [Notebook Outline](#toc)
 %% Cell type:markdown id:eb4bdea1-384f-41b3-8d8c-7bd568ae1537 tags:
 ## Additional Notes <a class="anchor" id="section_4"></a>
 %% Cell type:markdown id:cd25cb76-31eb-4c06-9bf8-a0407967b141 tags:
 1. Simulation data <br>
 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.
 2. Statistics <br>
 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.
 3. Scripts for visualizing network dynamics <br>
 The scripts for computing statistics and plotting the figures in S3 can be found in `./figures/MAM2EBRAINS`.
 %% Cell type:markdown id:ae940386-f9ec-4556-93bb-f63c7053cbd3 tags:
 Go back to [Notebook Outline](#toc)