diff --git a/.ipynb_checkpoints/multi-area-model-checkpoint.ipynb b/.ipynb_checkpoints/multi-area-model-checkpoint.ipynb
index 8ae4670dd7132a62d18c5cd80bb1bf6a702cca5f..e766a472f4359c459c13afe732002d489e9556b0 100644
--- a/.ipynb_checkpoints/multi-area-model-checkpoint.ipynb
+++ b/.ipynb_checkpoints/multi-area-model-checkpoint.ipynb
@@ -27,11 +27,10 @@
     "    * [2.2. Predict firing rates from theory](#section_2_2)\n",
     "    * [2.3. Extract and visualize interareal connectivity](#section_2_3)\n",
     "    * [2.4. Run a simulation](#section_2_4)\n",
-    "* [S3. Data Loading and Processing](#section_3)\n",
-    "* [S4. Simulation Results Visualization](#section_4) \n",
-    "    * [4.1. Instantaneous and mean firing rate across all populations](#section_4_1)\n",
-    "    * [4.2 Resting state plots](#section_4_2)\n",
-    "    * [4.3 Time-averaged population rates](#section_4_3)"
+    "* [S3. Simulation Results Visualization](#section_3) \n",
+    "    * [3.1. Instantaneous and mean firing rate across all populations](#section_3_1)\n",
+    "    * [3.2. Resting state plots](#section_3_2)\n",
+    "    * [3.3. Time-averaged population rates](#section_3_3)"
    ]
   },
   {
@@ -46,7 +45,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 25,
+   "execution_count": 49,
    "id": "9d6cc7d9-3110-4d96-9f9a-9ec7dee6d145",
    "metadata": {},
    "outputs": [],
@@ -64,7 +63,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 27,
+   "execution_count": 53,
    "id": "96517739",
    "metadata": {
     "tags": []
@@ -77,47 +76,30 @@
     "import nest\n",
     "import json\n",
     "import sys\n",
-    "from io import StringIO\n",
+    "from IPython.display import display, HTML\n",
     "\n",
     "from multiarea_model import MultiAreaModel\n",
     "from config import base_path, data_path\n",
     "\n",
-    "sys.path.append('./figures/MAM2EBRAINS')\n",
-    "\n",
-    "# Ignore and don't display warning messages\n",
-    "import warnings\n",
-    "warnings.filterwarnings(\"ignore\")\n",
-    "\n",
-    "# Redirect stdout to a dummy stream\n",
-    "original_stdout = sys.stdout\n",
-    "sys.stdout = StringIO()"
+    "sys.path.append('./figures/MAM2EBRAINS')"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 54,
    "id": "7e07b0d0",
    "metadata": {
     "tags": []
    },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Requirement already satisfied: nested_dict in /srv/main-spack-instance-2305/spack/var/spack/environments/ebrains-23-06/.spack-env/._view/ugqhkplzbs5vx62qrzwub7ficqgepl3r/lib/python3.8/site-packages (1.61)\n",
-      "Requirement already satisfied: dicthash in /srv/main-spack-instance-2305/spack/var/spack/environments/ebrains-23-06/.spack-env/._view/ugqhkplzbs5vx62qrzwub7ficqgepl3r/lib/python3.8/site-packages (0.0.2)\n",
-      "Requirement already satisfied: future in /srv/main-spack-instance-2305/spack/var/spack/environments/ebrains-23-06/.spack-env/._view/ugqhkplzbs5vx62qrzwub7ficqgepl3r/lib/python3.8/site-packages (from dicthash) (0.18.2)\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
-    "!pip install nested_dict dicthash;"
+    "%%capture captured\n",
+    "!pip install nested_dict dicthash"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 55,
    "id": "1d440c07-9b69-4e52-8573-26b13493bc5a",
    "metadata": {
     "tags": []
@@ -141,14 +123,12 @@
    ],
    "source": [
     "# Jupyter notebook display format setting\n",
-    "from IPython.display import display, HTML\n",
     "style = \"\"\"\n",
     "<style>\n",
     "table {float:left}\n",
     "</style>\n",
     "\"\"\"\n",
-    "display(HTML(style))\n",
-    "\n"
+    "display(HTML(style))"
    ]
   },
   {
@@ -209,7 +189,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 56,
    "id": "60265d52",
    "metadata": {},
    "outputs": [],
@@ -246,7 +226,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 18,
+   "execution_count": 57,
    "id": "6e4bed8d",
    "metadata": {},
    "outputs": [],
@@ -329,18 +309,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 19,
+   "execution_count": 58,
    "id": "ab25f9f8",
    "metadata": {},
    "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Initializing network from dictionary.\n",
-      "RAND_DATA_LABEL 3631\n"
-     ]
-    },
     {
      "name": "stderr",
      "output_type": "stream",
@@ -348,37 +320,10 @@
       "Error in library(\"aod\") : there is no package called ‘aod’\n",
       "Execution halted\n"
      ]
-    },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "No R installation or IndexError, taking hard-coded SLN fit parameters.\n",
-      "\n",
-      "\n",
-      "========================================\n",
-      "Customized parameters\n",
-      "--------------------\n",
-      "{'K_scaling': 0.005,\n",
-      " 'N_scaling': 0.005,\n",
-      " 'connection_params': {'K_stable': 'K_stable.npy',\n",
-      "                       'av_indegree_V1': 3950.0,\n",
-      "                       'fac_nu_ext_5E': 1.125,\n",
-      "                       'fac_nu_ext_6E': 1.41666667,\n",
-      "                       'fac_nu_ext_TH': 1.2,\n",
-      "                       'g': -11.0,\n",
-      "                       'replace_non_simulated_areas': 'het_poisson_stat'},\n",
-      " 'fullscale_rates': 'tests/fullscale_rates.json',\n",
-      " 'input_params': {'rate_ext': 10.0},\n",
-      " 'neuron_params': {'V0_mean': -150.0, 'V0_sd': 50.0}}\n",
-      "========================================\n",
-      "Simulation label: 27d81076e6d6e9e591684be053078477\n",
-      "Copied files.\n",
-      "Initialized simulation class.\n"
-     ]
     }
    ],
    "source": [
+    "%%capture captured\n",
     "M = MultiAreaModel(network_params, \n",
     "                   simulation=True,\n",
     "                   sim_spec=sim_params,\n",
@@ -396,19 +341,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 20,
+   "execution_count": 59,
    "id": "6a7ddf0e",
    "metadata": {},
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Iteration: 0\n",
-      "Mean-field theory predicts an average firing rate of 29.588 spikes/s across all populations.\n"
-     ]
-    }
-   ],
+   "outputs": [],
    "source": [
     "p, r = M.theory.integrate_siegert()\n",
     "print(\"Mean-field theory predicts an average \"\n",
@@ -433,7 +369,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 21,
+   "execution_count": 60,
    "id": "6316ac24",
    "metadata": {},
    "outputs": [],
@@ -449,7 +385,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 22,
+   "execution_count": 61,
    "id": "8408d463-557b-481b-afc1-5fbbbd67306d",
    "metadata": {},
    "outputs": [],
@@ -466,7 +402,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 23,
+   "execution_count": 62,
    "id": "445a722a",
    "metadata": {},
    "outputs": [],
@@ -483,18 +419,10 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 24,
+   "execution_count": 63,
    "id": "05512922-26e5-425f-90a4-0df7c2279ccf",
    "metadata": {},
    "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "Initializing network from dictionary.\n",
-      "RAND_DATA_LABEL 9405\n"
-     ]
-    },
     {
      "name": "stderr",
      "output_type": "stream",
@@ -503,35 +431,6 @@
       "Execution halted\n"
      ]
     },
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "No R installation or IndexError, taking hard-coded SLN fit parameters.\n",
-      "\n",
-      "\n",
-      "========================================\n",
-      "Customized parameters\n",
-      "--------------------\n",
-      "{'K_scaling': 1,\n",
-      " 'N_scaling': 1,\n",
-      " 'connection_params': {'K_stable': 'K_stable.npy',\n",
-      "                       'av_indegree_V1': 3950.0,\n",
-      "                       'fac_nu_ext_5E': 1.125,\n",
-      "                       'fac_nu_ext_6E': 1.41666667,\n",
-      "                       'fac_nu_ext_TH': 1.2,\n",
-      "                       'g': -11.0,\n",
-      "                       'replace_non_simulated_areas': 'het_poisson_stat'},\n",
-      " 'fullscale_rates': 'tests/fullscale_rates.json',\n",
-      " 'input_params': {'rate_ext': 10.0},\n",
-      " 'neuron_params': {'V0_mean': -150.0, 'V0_sd': 50.0}}\n",
-      "========================================\n",
-      "Simulation label: 155470013b00dadc9c4a4af26ef5090e\n",
-      "Copied files.\n",
-      "Initialized simulation class.\n",
-      "6.8556 6.052848304676828\n"
-     ]
-    },
     {
      "ename": "ValueError",
      "evalue": "The number of FixedLocator locations (33), usually from a call to set_ticks, does not match the number of ticklabels (32).",
@@ -539,30 +438,19 @@
      "traceback": [
       "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
       "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
-      "Cell \u001b[0;32mIn [24], line 2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mM2E_visualize_interareal_connectivity\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m visualize_interareal_connectivity\n\u001b[0;32m----> 2\u001b[0m visualize_interareal_connectivity(M)\n",
+      "Cell \u001b[0;32mIn [63], line 2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mM2E_visualize_interareal_connectivity\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m visualize_interareal_connectivity\n\u001b[0;32m----> 2\u001b[0m visualize_interareal_connectivity(M)\n",
       "File \u001b[0;32m~/MAM2EBRAINS/./figures/MAM2EBRAINS/M2E_visualize_interareal_connectivity.py:236\u001b[0m, in \u001b[0;36mvisualize_interareal_connectivity\u001b[0;34m(M)\u001b[0m\n\u001b[1;32m    233\u001b[0m X, Y \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mmeshgrid(x, y)\n\u001b[1;32m    235\u001b[0m ax\u001b[38;5;241m.\u001b[39mset_xticks([i \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m0.5\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m np\u001b[38;5;241m.\u001b[39marange(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;28mlen\u001b[39m(area_list) \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m1\u001b[39m)])\n\u001b[0;32m--> 236\u001b[0m \u001b[43max\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mset_xticklabels\u001b[49m\u001b[43m(\u001b[49m\u001b[43marea_list\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrotation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m90\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msize\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m6.\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m    238\u001b[0m ax\u001b[38;5;241m.\u001b[39mset_yticks([i \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m0.5\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m np\u001b[38;5;241m.\u001b[39marange(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;28mlen\u001b[39m(area_list) \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m1\u001b[39m)])\n\u001b[1;32m    239\u001b[0m ax\u001b[38;5;241m.\u001b[39mset_yticklabels(area_list[::\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m], size\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m6.\u001b[39m)\n",
       "File \u001b[0;32m/srv/main-spack-instance-2305/spack/opt/spack/linux-ubuntu20.04-x86_64/gcc-10.3.0/py-matplotlib-3.6.2-lhkot3cmeebfk5dp74dnubweq56upksc/lib/python3.8/site-packages/matplotlib/axes/_base.py:73\u001b[0m, in \u001b[0;36m_axis_method_wrapper.__set_name__.<locals>.wrapper\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m     72\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mwrapper\u001b[39m(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m---> 73\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mget_method\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n",
       "File \u001b[0;32m/srv/main-spack-instance-2305/spack/opt/spack/linux-ubuntu20.04-x86_64/gcc-10.3.0/py-matplotlib-3.6.2-lhkot3cmeebfk5dp74dnubweq56upksc/lib/python3.8/site-packages/matplotlib/axis.py:1968\u001b[0m, in \u001b[0;36mAxis._set_ticklabels\u001b[0;34m(self, labels, fontdict, minor, **kwargs)\u001b[0m\n\u001b[1;32m   1966\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m fontdict \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m   1967\u001b[0m     kwargs\u001b[38;5;241m.\u001b[39mupdate(fontdict)\n\u001b[0;32m-> 1968\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mset_ticklabels\u001b[49m\u001b[43m(\u001b[49m\u001b[43mlabels\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mminor\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mminor\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n",
       "File \u001b[0;32m/srv/main-spack-instance-2305/spack/opt/spack/linux-ubuntu20.04-x86_64/gcc-10.3.0/py-matplotlib-3.6.2-lhkot3cmeebfk5dp74dnubweq56upksc/lib/python3.8/site-packages/matplotlib/axis.py:1890\u001b[0m, in \u001b[0;36mAxis.set_ticklabels\u001b[0;34m(self, ticklabels, minor, **kwargs)\u001b[0m\n\u001b[1;32m   1886\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28misinstance\u001b[39m(locator, mticker\u001b[38;5;241m.\u001b[39mFixedLocator):\n\u001b[1;32m   1887\u001b[0m     \u001b[38;5;66;03m# Passing [] as a list of ticklabels is often used as a way to\u001b[39;00m\n\u001b[1;32m   1888\u001b[0m     \u001b[38;5;66;03m# remove all tick labels, so only error for > 0 ticklabels\u001b[39;00m\n\u001b[1;32m   1889\u001b[0m     \u001b[38;5;28;01mif\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(locator\u001b[38;5;241m.\u001b[39mlocs) \u001b[38;5;241m!=\u001b[39m \u001b[38;5;28mlen\u001b[39m(ticklabels) \u001b[38;5;129;01mand\u001b[39;00m \u001b[38;5;28mlen\u001b[39m(ticklabels) \u001b[38;5;241m!=\u001b[39m \u001b[38;5;241m0\u001b[39m:\n\u001b[0;32m-> 1890\u001b[0m         \u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\n\u001b[1;32m   1891\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mThe number of FixedLocator locations\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m   1892\u001b[0m             \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m (\u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mlen\u001b[39m(locator\u001b[38;5;241m.\u001b[39mlocs)\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m), usually from a call to\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m   1893\u001b[0m             \u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m set_ticks, does not match\u001b[39m\u001b[38;5;124m\"\u001b[39m\n\u001b[1;32m   1894\u001b[0m             \u001b[38;5;124mf\u001b[39m\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124m the number of ticklabels (\u001b[39m\u001b[38;5;132;01m{\u001b[39;00m\u001b[38;5;28mlen\u001b[39m(ticklabels)\u001b[38;5;132;01m}\u001b[39;00m\u001b[38;5;124m).\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m   1895\u001b[0m     tickd \u001b[38;5;241m=\u001b[39m {loc: lab \u001b[38;5;28;01mfor\u001b[39;00m loc, lab \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mzip\u001b[39m(locator\u001b[38;5;241m.\u001b[39mlocs, ticklabels)}\n\u001b[1;32m   1896\u001b[0m     func \u001b[38;5;241m=\u001b[39m functools\u001b[38;5;241m.\u001b[39mpartial(\u001b[38;5;28mself\u001b[39m\u001b[38;5;241m.\u001b[39m_format_with_dict, tickd)\n",
       "\u001b[0;31mValueError\u001b[0m: The number of FixedLocator locations (33), usually from a call to set_ticks, does not match the number of ticklabels (32)."
      ]
-    },
-    {
-     "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfoAAAGqCAYAAADjm4a+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAAAe+0lEQVR4nO3de7RedX3n8feXhAhyEWvCPTHRQTF1ocIxMlQplykmGVtqx85wGR0Yl5EqjrPWXKROl1odR+20S2vFxmgz0JnBTKuMoKJYrIBdgCXhftVIgISIIdyUa0jynT9+v8dsHw/kAOfsc85v3q+1nkX29ffdez/n+ezf3vt5iMxEkiS1aZfJLkCSJE0cg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWrYToM+IuZHRNbX9oi4LyLOi4g9+yhQkiQ9d8+mR38tcCpwFXAy8AcTUpEkSRo3zybo7wMuAdbU4ZnjX44kSRpPsbNfxouI+cC6odEbgddl5qYJqkuSJI2DZ9Oj/wHwW8DHgQOB90xIRZIkadw8m6DfnJmXUIIeYOkE1CNJksbRs7nPfmBEnAS8qQ7fOf7lSJKk8fRc7tH/DLgSeHdm3jVxpUmSpOdrp0EvSZKmL38ZT5Kkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkN6zXoI2JlRGyKiJueZnpExGcjYm1E3BARh/dZnyRJrem7R38OsPgZpi8BDqmvZcBf9lCTJEnN6jXoM/Ny4IFnmOVE4K+zuArYJyIO6Kc6SZLaM3OyCxhyELC+M7yhjvvJ8IwRsYzS62ePPfY44tBDD+2lQKlFa9as2ZyZcya7Dknjb6oFfYwyLkebMTNXACsARkZGcvXq1RNZl9S0iLhrsmuQNDGm2lP3G4C5neGDgY2TVIskSdPeVAv6C4F31KfvjwQezsxfuWwvSZLGptdL9xHxZeAYYHZEbAA+DOwKkJnLgYuApcBa4DHg9D7rkySpNb0GfWaevJPpCby3p3IkSWreVLt0L0mSxpFBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqmEEvSVLDDHpJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqmEEvSVLDDHpJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqWO9BHxGLI+L2iFgbEWeNMv1FEfH1iLg+Im6OiNP7rlGSpFb0GvQRMQM4G1gCLAROjoiFQ7O9F7glM18DHAP8WUTM6rNOSZJa0XePfhGwNjPvyMwtwCrgxKF5EtgrIgLYE3gA2NpvmZIktaHvoD8IWN8Z3lDHdX0OeBWwEbgReH9mbh9eUUQsi4jVEbH6vvvum6h6JUma1voO+hhlXA4Nvxm4DjgQeC3wuYjY+1cWylyRmSOZOTJnzpzxrlOSpCb0HfQbgLmd4YMpPfeu04Hzs1gLrAMO7ak+SZKa0nfQXw0cEhEL6gN2JwEXDs1zN3A8QETsB7wSuKPXKiVJasTMPhvLzK0RcSZwMTADWJmZN0fEGXX6cuBjwDkRcSPlUv8HMnNzn3VKktSKXoMeIDMvAi4aGre88++NwAl91yVJUov8ZTxJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqmEEvSVLDDHpJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqmEEvSVLDDHpJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGtZ70EfE4oi4PSLWRsRZTzPPMRFxXUTcHBGX9V2jJEmtmNlnYxExAzgb+C1gA3B1RFyYmbd05tkH+DywODPvjoh9+6xRkqSW9N2jXwSszcw7MnMLsAo4cWieU4DzM/NugMzc1HONkiQ1o++gPwhY3xneUMd1vQJ4cURcGhFrIuIdvVUnSVJjer10D8Qo43JoeCZwBHA8sDtwZURclZk//KUVRSwDlgHMmzdvAkqVJGn667tHvwGY2xk+GNg4yjzfzsxHM3MzcDnwmuEVZeaKzBzJzJE5c+ZMWMGSJE1nfQf91cAhEbEgImYBJwEXDs1zAfCmiJgZES8E3gDc2nOdkiQ1oddL95m5NSLOBC4GZgArM/PmiDijTl+embdGxLeBG4DtwJcy86Y+65QkqRWROXyLfPoZGRnJ1atXT3YZ0rQVEWsyc2Sy65A0/vxlPEmSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqmEEvSVLDDHpJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9NI0dOeddxIRRAS77bYbc+fO5dRTT2XdunWTXZqkKcagl6ax173udSxfvpzjjjuO8847j6OOOopNmzZNdlmSphCDXprGDjzwQE477TTOPfdc3vWud3HvvffyhS98YbLLkjSFGPRSI5YsWQLA9ddfP8mVSJpKDHqpEZkJQERMciWSphKDXmrExRdfDMBhhx02yZVImkpmTnYBkp67jRs3cs4553DZZZdx7rnnsv/++7Ns2bIJaSsiVgJvATZl5qsnpBFJ484evTSNXXvttbz73e/mkksu4ZRTTuHKK69kv/32m6jmzgEWT9TKJU0Me/TSNDR//vxf3JPvS2ZeHhHze21U0vNm0EuaELNnz8758+dPdhnStLVmzZrNmTnn+a7HoJc0IebPn8/q1asnuwxp2oqIu8ZjPd6jlySpYQa9pDGJiMXAd4FDIuKsya5H0tgY9JJ2KiJmAH8L7EX53PhoRPzR5FYlaSwMekljsQi4IjP3zcxdgQ8D2ya5JkljYNBLGouDgPWd4Q11nKQpzqfuJY3FaD+g/4xf5F+zZs3m8XpqWPr/1EvHYyUGvaSx2ADM7QwfDGx8pgXG4/u/kp4/L91LGourKU/bL4iIWcBJwIWTXJOkMbBHL2mnMnNrRJwJXAzMAFZm5s2TXJakMTDoJY1JZl4EXDTZdUh6drx0L0lSwwx6SZIa5qV7SeOm/kzun1Pu438pMz9Zx68E3gI8CvwasK4ucglwBPAyypP8PwE2Ux70Wwy8ANgbOAC4rS5zfmZ+NCKOA/4OeKSu73zg3wI/B3YH5gM3DZYBPgtcALwReAq4G/gW8Jt1nj2Bl1O+TXBfrWUe8ELKbwbcSvlK4fnAw8C/r9OfAGYBD2TmfnV7f7u2FcB24FOZ+cGI+H3gI8BC4MHa5u2UHx/aH9itjtultvcUsG8dv1dt66H62gZsret/fZ13sK59gRfVumbUbX2ozn8u8F+BPWp9t9Rl5tXt25uSDU/U9e1a9+khdRr1WGyr054AXln305Z6/O4FTgMur9szE/hpfe1V63gFcGfdt5uBTcCBta4ZtfbBMvt31rO1buuGWkPUZV5exw/WtT/w4rpNs+r+Xl+XvwD4w9rOoLZNlPfMjLpdL6jru6UOP1n3wcz6uq8uM4vyjZSZdZnN7PhGyoJ67J6qtd5dt2NwrH9U23yQ8n7bvx6nLXV77qzLbs3MkbpOIuLOeky2DU8bVWZO+9cRRxyRkp47YHU+z7/D+gH5Y0pozwKuBxbWaUcDh1MC+RudZQ6o44+hhO4PKSEYwJ51nuPrB+GRQ+2dDdwztL47gdl1fd8Ymv9c4L8D36j17TM0/ThKaL2UEj7rKCcMx9QP7tPqfK+mnEC8GfgmsBo4Ebips64v11oOp4Tep+r4V1FC8Vrg1PqBPrtOOwE4ti6zDfhsHb93Z/9trB/wg2XmAv9ICYnuuj5S98/hQ+OPpZxcHT88beg4bQM+Wcd9B1hSp72PcmIxWN/VwFvrMpspJ3mDY3h2ZxserK+FdR+8se63EeAuYG2ddhLw+rrMw51lDgEOr+P/U63hqDp8BHBFXc/6zrr+FPh0nafbxrGUE5A31GndZQ7otPMw8EAd/x3glLqdSynvk8EyVwOL6zKb67E4EvgT4EN1/AcpJxhvZsd74Pt1++8E1tRlTgBmdtr/7NP8rd3ZPW47e3npXtJ4WQSszcw7MnMLsIoSgGTm5ZQPzV+SmT/JzGvq4DZKL/agev7xSB0/k9LD/cUP9ETEwZQPxrvHUlhEDMLym7XdLZn50NBshwOPZebgR35mUoJ+Rn0NemmvAq6i9PCS0js8cmhdv0HprT1A6aH9bm331sy8nfIh/rOhffGdzPxeXSYpvVsy82ed/fdL+wH4NPAuRv854nX86j7/A0qAf3eUaQy1c/5gNLB3nTYYHngl8LV6DB+nXIW5lXKi9M+AT9T5HqH0ig+q++Af6rjB+n5Up63KzKvr+CfrOg/KzB913icz6/YeUIf/sG5XdtdV13/PcBt13o9m5g/qtG2d9n+SmddExOAqwY11mQS21RpeRAntwfpeSfk2CpQTgMFVhBOBL9bxqyjvoey8BwbHLOo2ZX0PbO1s/4GMh+d7Fj8VXvbopeeH8enRv41yuX4w/Hbgc53h+ZTwuZ/S2/8W8Ot12jGU3tsWSu/p1+sH43WUD/vHu8sAXwGWsSMwB+PXAddQLjk/0lnm9yg9329TPqTvH7TTqe8iyq2FwTKfqOt4qH7oDsb/NqXX+jt1PY9SLiP/qLOuR+q0Wyg98IeH9tWldT1bar1rgGWd/bS9bsuaup0fp5xoPEHpnV5D6VF+ty6zYWhdX6f0+m6lBMr1dfx64I+BH1BOVp4apf1/WdsfjP8Q5YRqPeXqRLedHwMn1uXur21tAf4d8FBne9dTwu/aTjuXUnq06+syv5hWpz9ej9W1nX2wvh7brZT3Rncf3NnZD9d29sEN9XgMxnf3wWWdberWdnQ95oPx3X1wT30Npv2YclXjurqN2+u+eYwd7+FHhvbpsrr9t9fx93aPQa3hMTrvgaH3zzqGjtsz/m0+3z/uqfAy6KXnh/EJ+t/nV4P+LzrD8ynBN7gkv5QajpTe2bWUQP7F+DptLuVS66vrtHuAz1MuaS+hXIpfSulhHViXeRmlN3Z0nXZ3DYdjKfdI/5xyeX3Q/ixKUC3IHWH3GDCH8kzB14F/3WnnnfUD/B+A5ZQe3ZOdmh+q7cynhM22oX11KSXob6vD+1LC+GjK7YWfU3p63fHzayD8GeV++JrONv7iRKMuczPl5GkB5XLyyjr+ccpJUlB6nFtGaed/AZs667oP+KM6/F7g0c602ygnUGvqPt1KuTd/PfDzOt+e7AjbbjuXAm+q035vaNofU07ihsfvSTn+F4yyD+6q2z1YZrAP9qLcg//uKPvgN+s+GG7ni5QTq8H47j54ez0+g2ndffAXlJOnoyknPUfXZfap41/daedayonOocD36jKD9v8L5aTyl45N5/1z4PD75pn+Nr10L2m8jOVncrdnvSSf5Xv5u0bE/sA5wP/MzPM742fX+dZTPqQX12l7UC6F3wD8D8q99VMol4a31GXuAL4KLKrLBLAxM79X2/8K5URh0M4SYE1mrqt1JuWDOjPzAcr/oveoQW3ABZn52sx8I+VS9/eBGNTMjofOoPTCt3emdW2t9W4C/i/lsvLxwN31HGwwflGd/2HKfd6X1309F/ga5SGul0bE/nWZv6WESFJOYBbV8fcA99WTu+vr9Nmddo6s63+gU9fgIUAotz5260xbBfxNXe7wuo/Pqet6PCLm1uPwVeCeoe0Jyn3sc+pxH0x7D/B+4L8NjT+yrmdl3f7hfTCXcqn7is4+WEQ51iuA/br7gHK5/IOUK0Lf77RzFPAO4Iud9vcCnoiIXSknfLt2pq0C/iYzj6AE/U9ruw/WYwnlFtDDlPfwoJ296368jXLSs6jzHngL8C+e5j1AZm7sHINfmjYag17SeBnLz+TOrPc/iYhFlM+gP6H0xj7dGT94upqImE+533tbnfYzyr3RI2sbf095on4G5XIrEbGA8mDTTXWZBO6KiN+o7R/PjnvR9wMnA98c1EbpOb4AeKyeiBwP3NqpeZeI2D8i5lF6dnfU5e6v/70E+Df133MovcjBtIEX1HUREXtQriIcRXngLTrjf4cd3x7Ym/IcxI2U3votlJOce4F1mXlvXeaf12V2p9xTvqmOD3Y8Of+q2v7mOu0EyknMOsrl5EH729mRFcdST6bqtKW1nXPrdn6ks67LKD3vH9X9f0Fn2k2Ue9vrqfexO/vgLZRA//zQPjiBcivioVH2wTWUKwavHdoHJ9T2Hx1lH/wVJfCfGNoH/4pyFeRPR9kHf1WP5a1Dta2LiD2Bj7HjkvtTlJMRgDMo783bOvvtcWC3iJhDeX+v67wHTqK8nxnaZ0TEHhGx12jTnk7U7v+0NjIykqtXr57sMqRpKyLW5M6+ojO29SwFPsOOn8n9eEScQbnUfTDlUuPgnuRTwJWUoNpIuUQOJfSvpITrQ5Sn6GdSnmbeC/hKZv7n+pO8/wF4CeXD/nrKh/v9tZ1dajt7UXp151GCZz/Kh+4PKZdM11OeUv8EcDrlBGOv2t5LKKHwQnZ8NeorwD8F/gnlwaunKKE9g3IS8h/rch9nx1eYNwNnUQLvbMqtgqivpygnHb9W/z2rLjfYT49QbkVQ29hWxz9Wx11DCYhdOuvavba1Z11mK+WrYP+b8rT/HvU1s7PMw5Sgn1OXe4oShN+j3GbYzo6vuA2Wibov59Y2Bg+SfZ0SwiN13hmUE4RNlKft31K3c/D1uC31uO1X1/9knb6d0gOP2sYTtY3B++Gpuq1z2fEVNuq69qnb+CQ7rvZsAv4P5TjPrtMGX4G8v/53LqVX/qK6rk2UKza/W9f3RF1mcLVk9zp+l/raXMdfSwnswdcgH6nHbTbl/Tmjthd1Pz9c99lPKbcGFlBOUH5a9+dhmbk0Il5G6cVT98N5mflxnoFBL2ncgl7S1OOle0mSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqmEEvSVLDDHpJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ3rPegjYnFE3B4RayPirGeY7/URsS0i3tZnfZIktaTXoI+IGcDZwBJgIXByRCx8mvk+BVzcZ32SJLWm7x79ImBtZt6RmVuAVcCJo8z3PuCrwKY+i5MkqTV9B/1BwPrO8IY67hci4iDgrcDyHuuSJKlJfQd9jDIuh4Y/A3wgM7c944oilkXE6ohYfd99941XfZIkNWVmz+1tAOZ2hg8GNg7NMwKsigiA2cDSiNiamV/rzpSZK4AVACMjI8MnC5Ikif6D/mrgkIhYANwDnASc0p0hMxcM/h0R5wDfGA55SZI0Nr0GfWZujYgzKU/TzwBWZubNEXFGne59eUmSxlHfPXoy8yLgoqFxowZ8Zp7WR02SJLXKX8aTJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqmEEvSVLDDHpJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqmEEvSVLDDHpJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWG9B31ELI6I2yNibUScNcr0UyPihvq6IiJe03eNkiS1otegj4gZwNnAEmAhcHJELByabR3wm5l5GPAxYEWfNUqS1JK+e/SLgLWZeUdmbgFWASd2Z8jMKzLzwTp4FXBwzzVKktSMvoP+IGB9Z3hDHfd03gl8a0IrkiSpYTN7bi9GGZejzhhxLCXo3/g005cBywDmzZs3XvVJktSUvnv0G4C5neGDgY3DM0XEYcCXgBMz8/7RVpSZKzJzJDNH5syZMyHFSpI03fUd9FcDh0TEgoiYBZwEXNidISLmAecDb8/MH/ZcnyRJTen10n1mbo2IM4GLgRnAysy8OSLOqNOXAx8CXgJ8PiIAtmbmSJ91SpLUisgc9Rb5tDIyMpKrV6+e7DKkaSsi1nhCLbXJX8aTJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqmEEvSVLDDHpJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWEGvSRJDTPoJUlqmEEvSVLDDHpJkhpm0EuS1DCDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIaZtBLktQwg16SpIYZ9JIkNcyglySpYQa9JEkNM+glSWqYQS9JUsMMekmSGmbQS5LUMINekqSGGfSSJDXMoJckqWG9B31ELI6I2yNibUScNcr0iIjP1uk3RMThfdcoSVIreg36iJgBnA0sARYCJ0fEwqHZlgCH1Ncy4C/7rFGSpJb03aNfBKzNzDsycwuwCjhxaJ4Tgb/O4ipgn4g4oOc6JUlqQt9BfxCwvjO8oY57tvNIkqQxmNlzezHKuHwO8xARyyiX9gGejIibnmdtfZgNbJ7sIsbAOsfPdKgR4JWTXYCkidF30G8A5naGDwY2Pod5yMwVwAqAiFidmSPjW+r4s87xNR3qnA41QqlzsmuQNDH6vnR/NXBIRCyIiFnAScCFQ/NcCLyjPn1/JPBwZv6k5zolSWpCrz36zNwaEWcCFwMzgJWZeXNEnFGnLwcuApYCa4HHgNP7rFGSpJb0femezLyIEubdccs7/07gvc9ytSvGobQ+WOf4mg51TocaYfrUKelZipKrkiSpRf4EriRJDZtWQT9dfj53DHWeWuu7ISKuiIjXTLUaO/O9PiK2RcTb+qyv0/5O64yIYyLiuoi4OSIu67vGWsPOjvmLIuLrEXF9rbP3Z08iYmVEbHq6r6JOlb8fSeMsM6fFi/Lw3o+BlwGzgOuBhUPzLAW+Rfku/pHAD6ZonUcBL67/XtJ3nWOpsTPf31OeqXjbFN2X+wC3APPq8L5TtM4PAp+q/54DPADM6rnOo4HDgZueZvqk//348uVr/F/TqUc/XX4+d6d1ZuYVmflgHbyK8lsBU6rG6n3AV4FNfRbXMZY6TwHOz8y7ATJzMmodS50J7BURAexJCfqtfRaZmZfXdp/OVPj7kTTOplPQT5efz322NbyT0ovq005rjIiDgLcCy5k8Y9mXrwBeHBGXRsSaiHhHb9XtMJY6Pwe8ivLjTzcC78/M7f2UN2ZT4e9H0jjr/et1z8O4/XzuBBtzDRFxLCXo3zihFY3S9Cjjhmv8DPCBzNxWOqGTYix1zgSOAI4HdgeujIirMvOHE11cx1jqfDNwHXAc8HLg7yLi+5n5swmu7dmYCn8/ksbZdAr6cfv53Ak2phoi4jDgS8CSzLy/p9oGxlLjCLCqhvxsYGlEbM3Mr/VSYTHWY745Mx8FHo2Iy4HXAH0G/VjqPB34ZGYmsDYi1gGHAv/YT4ljMhX+fiSNs+l06X66/HzuTuuMiHnA+cDbe+55jrnGzFyQmfMzcz7wFeA9PYf8mOoELgDeFBEzI+KFwBuAW6dgnXdTrjoQEftR/icyd/Ra5c5Nhb8fSeNs2vToc5r8fO4Y6/wQ8BLg87XHvDV7/B+fjLHGSTeWOjPz1oj4NnADsB34Umb2+n8yHOP+/BhwTkTcSLlE/oHM7PX/ahcRXwaOAWZHxAbgw8CunRon/e9H0vjzl/EkSWrYdLp0L0mSniWDXpKkhhn0kiQ1zKCXJKlhBr0kSQ0z6CVJaphBL0lSwwx6SZIa9v8Anz2BVULhdqAAAAAASUVORK5CYII=\n",
-      "text/plain": [
-       "<Figure size 493.603x435.805 with 2 Axes>"
-      ]
-     },
-     "metadata": {
-      "needs_background": "light"
-     },
-     "output_type": "display_data"
     }
    ],
    "source": [
-    "from M2E_visualize_interareal_connectivity import visualize_interareal_connectivity\n",
-    "visualize_interareal_connectivity(M);"
+    "%%capture captured\n",
+    "# from M2E_visualize_interareal_connectivity import visualize_interareal_connectivity\n",
+    "# visualize_interareal_connectivity(M)"
    ]
   },
   {
@@ -687,67 +575,6 @@
     "Go back to [Notebook structure](#toc)"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "id": "57ff902c-d6ce-4f96-9e4f-8e3e7166ab66",
-   "metadata": {
-    "tags": []
-   },
-   "source": [
-    "## S3. Data Loading and Processing <a class=\"anchor\" id=\"section_3\"></a>"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 12,
-   "id": "f5b58845-4d1a-430f-83f4-402fdf918aef",
-   "metadata": {
-    "tags": []
-   },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "loading spikes\n",
-      "Loading data from file\n",
-      "Computing population rates done\n",
-      "Loading data from file\n",
-      "Computing synchrony done\n",
-      "Loading data from file\n",
-      "Computing population LvR done\n",
-      "Loading data from file\n",
-      "Loading data from file\n",
-      "Computing rate time series done\n",
-      "pop_LvR\n",
-      "pop_rates\n",
-      "synchrony\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "python3: can't open file './../Schmidt2018_dyn/compute_bold_signal.py': [Errno 2] No such file or directory\n"
-     ]
-    }
-   ],
-   "source": [
-    "label_spikes = M.simulation.label\n",
-    "label = M.simulation.label\n",
-    "\n",
-    "from MAM2EBRAINS_LOAD_DATA import load_data\n",
-    "A, tsteps, firing_rate = load_data(M)"
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "id": "2da9728d-4481-4a15-b810-d125e39cbe4e",
-   "metadata": {},
-   "source": [
-    "Go back to [Notebook structure](#toc)"
-   ]
-  },
   {
    "cell_type": "markdown",
    "id": "bb71c922",
@@ -755,7 +582,7 @@
     "tags": []
    },
    "source": [
-    "## S4. Simulation Results Visualziation <a class=\"anchor\" id=\"section_4\"></a>"
+    "## S3. Simulation Results Visualziation <a class=\"anchor\" id=\"section_3\"></a>"
    ]
   },
   {
@@ -765,7 +592,7 @@
     "tags": []
    },
    "source": [
-    "### 4.1. Instantaneous and mean firing rate across all populations <a class=\"anchor\" id=\"section_4_1\"></a>"
+    "### 3.1. Instantaneous and mean firing rate across all populations <a class=\"anchor\" id=\"section_3_1\"></a>"
    ]
   },
   {
@@ -788,7 +615,7 @@
     }
    ],
    "source": [
-    "from M2E_VISUALIZATION import plot_instan_mean_firing_rate\n",
+    "from M2E_visualize_instantaneous_and_mean_firing_rates import plot_instan_mean_firing_rate\n",
     "plot_instan_mean_firing_rate(M)"
    ]
   },
@@ -797,7 +624,7 @@
    "id": "e91c436e-db94-4cd7-a531-29c032efeeae",
    "metadata": {},
    "source": [
-    "### 4.2 Resting state plots <a class=\"anchor\" id=\"section_4_2\"></a>"
+    "### 3.2 Resting state plots <a class=\"anchor\" id=\"section_3_2\"></a>"
    ]
   },
   {
@@ -830,8 +657,8 @@
     }
    ],
    "source": [
-    "from MAM2EBRAINS_VISUALIZATION import plot_resting_state\n",
-    "plot_resting_state(M, A, label_spikes, data_path, sim_params)"
+    "from M2E_visualize_resting_state import plot_resting_state\n",
+    "plot_resting_state(M, A, label_spikes, data_path)"
    ]
   },
   {
@@ -841,7 +668,7 @@
     "tags": []
    },
    "source": [
-    "### 4.3 Time-averaged population rates <a class=\"anchor\" id=\"section_4_3\"></a>\n",
+    "### 3.3 Time-averaged population rates <a class=\"anchor\" id=\"section_4_3\"></a>\n",
     "Plot overview over time-averaged population rates encoded in colors with areas along x-axis and populations along y-axis."
    ]
   },
@@ -904,11 +731,51 @@
     }
    ],
    "source": [
+    "%%capture captured\n",
     "# 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']\n",
     "# output = {'pdf', 'png', 'eps'}, optional\n",
     "A.show_rates()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "15bcb6a1-9117-4538-8582-7a89e06167da",
+   "metadata": {},
+   "source": [
+    "0 V1\n",
+    "1 V2\n",
+    "2 VP\n",
+    "3 V3\n",
+    "4 PIP\n",
+    "5 V3A\n",
+    "6 MT\n",
+    "7 V4t\n",
+    "8 V4\n",
+    "9 PO\n",
+    "10 VOT\n",
+    "11 DP\n",
+    "12 MIP\n",
+    "13 MDP\n",
+    "14 MSTd\n",
+    "15 VIP\n",
+    "16 LIP\n",
+    "17 PITv\n",
+    "18 PITd\n",
+    "19 AITv\n",
+    "20 MSTl\n",
+    "21 FST\n",
+    "22 CITv\n",
+    "23 CITd\n",
+    "24 7a\n",
+    "25 STPp\n",
+    "26 STPa\n",
+    "27 FEF\n",
+    "28 46\n",
+    "29 TF\n",
+    "30 TH\n",
+    "31 AITd"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "ef74ca3e-98dc-49c9-a4a0-2c640e29b1d9",
diff --git a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_VISUALIZATION-checkpoint.py b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_VISUALIZATION-checkpoint.py
index b22bb9fd68ab69f79d77b40ec4498ef4b5ccc1b5..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644
--- a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_VISUALIZATION-checkpoint.py
+++ b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_VISUALIZATION-checkpoint.py
@@ -1,452 +0,0 @@
-import json
-import numpy as np
-import os
-
-import sys
-sys.path.append('./figures/Schmidt2018_dyn')
-
-from helpers import original_data_path, population_labels
-from multiarea_model import MultiAreaModel
-from plotcolors import myred, myblue
-
-import matplotlib.pyplot as pl
-from matplotlib import gridspec
-# from matplotlib import rc_file
-# rc_file('plotstyle.rc')
-
-icolor = myred
-ecolor = myblue
-
-# label_spikes = M.simulation.label
-label = M.simulation.label
-
-from MAM2EBRAINS_LOAD_DATA import load_and_create_data
-A = load_and_create_data(M)
-
-def plot_instan_mean_firing_rate(M):
-    # load spike data and calculate instantaneous and mean firing rates
-    data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + "-spikes-1-0.dat", skiprows=3)
-    tsteps, spikecount = np.unique(data[:,1], return_counts=True)
-    firing_rate = spikecount / M.simulation.params['dt'] * 1e3 / np.sum(M.N_vec)
-    
-    # visualize calculate instantaneous and mean firing rates
-    ax = pl.subplot()
-    ax.plot(tsteps, rate)
-    ax.plot(tsteps, np.average(rate)*np.ones(len(tsteps)), label='mean')
-    ax.set_title('Instantaneous and mean firing rate across all populations')
-    ax.set_xlabel('time (ms)')
-    ax.set_ylabel('firing rate (spikes / s)')
-    ax.set_xlim(0, sim_params['t_sim'])
-    ax.set_ylim(0, 50)
-    ax.legend()
-
-def set_boxplot_props(d):
-    for i in range(len(d['boxes'])):
-        if i % 2 == 0:
-            d['boxes'][i].set_facecolor(icolor)
-            d['boxes'][i].set_color(icolor)
-        else:
-            d['boxes'][i].set_facecolor(ecolor)
-            d['boxes'][i].set_color(ecolor)
-    pl.setp(d['whiskers'], color='k')
-    pl.setp(d['fliers'], color='k', markerfacecolor='k', marker='+')
-    pl.setp(d['medians'], color='none')
-    pl.setp(d['caps'], color='k')
-    pl.setp(d['means'], marker='x', color='k',
-            markerfacecolor='k', markeredgecolor='k', markersize=3.)
-
-def plot_resting_state(M, A, label_spikes, data_path, sim_params):
-    t_sim = sim_params["t_sim"]
-    
-    """
-    Figure layout
-    """
-
-    nrows = 4
-    ncols = 4
-    # width = 7.0866
-    width = 10
-    panel_wh_ratio = 0.7 * (1. + np.sqrt(5)) / 2.  # golden ratio
-
-    height = width / panel_wh_ratio * float(nrows) / ncols
-    pl.rcParams['figure.figsize'] = (width, height)
-
-
-    fig = pl.figure()
-    axes = {}
-
-    gs1 = gridspec.GridSpec(1, 3)
-    gs1.update(left=0.06, right=0.72, top=0.95, wspace=0.4, bottom=0.35)
-    # axes['A'] = pl.subplot(gs1[:-1, :1])
-    # axes['B'] = pl.subplot(gs1[:-1, 1:2])
-    # axes['C'] = pl.subplot(gs1[:-1, 2:])
-    axes['A'] = pl.subplot(gs1[:1, :1])
-    axes['B'] = pl.subplot(gs1[:1, 1:2])
-    axes['C'] = pl.subplot(gs1[:1, 2:])
-
-    gs2 = gridspec.GridSpec(3, 1)
-    gs2.update(left=0.78, right=0.95, top=0.95, bottom=0.35)
-    axes['D'] = pl.subplot(gs2[:1, :1])
-    axes['E'] = pl.subplot(gs2[1:2, :1])
-    axes['F'] = pl.subplot(gs2[2:3, :1])
-
-
-    gs3 = gridspec.GridSpec(1, 1)
-    gs3.update(left=0.1, right=0.95, top=0.3, bottom=0.075)
-    axes['G'] = pl.subplot(gs3[:1, :1])
-
-    areas = ['V1', 'V2', 'FEF']
-
-    labels = ['A', 'B', 'C']
-    for area, label in zip(areas, labels):
-        label_pos = [-0.2, 1.01]
-        # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label + ': ' + area,
-        #         fontdict={'fontsize': 10, 'weight': 'bold',
-        #                   'horizontalalignment': 'left', 'verticalalignment':
-        #                   'bottom'}, transform=axes[label].transAxes)
-        pl.text(label_pos[0], label_pos[1], label + ': ' + area,
-                 fontdict={'fontsize': 10, 'weight': 'bold', 
-                           'horizontalalignment': 'left', 'verticalalignment': 
-                           'bottom'}, transform=axes[label].transAxes)
-
-    label = 'G'
-    label_pos = [-0.1, 0.92]
-    # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
-    #         fontdict={'fontsize': 10, 'weight': 'bold',
-    #                   'horizontalalignment': 'left', 'verticalalignment':
-    #                   'bottom'}, transform=axes[label].transAxes)
-    pl.text(label_pos[0], label_pos[1], label,
-             fontdict={'fontsize': 10, 'weight': 'bold', 
-                       'horizontalalignment': 'left', 'verticalalignment': 
-                       'bottom'}, transform=axes[label].transAxes)
-
-    labels = ['E', 'D', 'F']
-    for label in labels:
-        label_pos = [-0.2, 1.05]
-        # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
-        #         fontdict={'fontsize': 10, 'weight': 'bold',
-        #                   'horizontalalignment': 'left', 'verticalalignment':
-        #                   'bottom'}, transform=axes[label].transAxes)
-        pl.text(label_pos[0], label_pos[1], label,
-             fontdict={'fontsize': 10, 'weight': 'bold', 
-                       'horizontalalignment': 'left', 'verticalalignment': 
-                       'bottom'}, transform=axes[label].transAxes)
-        
-
-    labels = ['A', 'B', 'C', 'D', 'E', 'F']
-
-    for label in labels:
-        axes[label].spines['right'].set_color('none')
-        axes[label].spines['top'].set_color('none')
-        axes[label].yaxis.set_ticks_position("left")
-        axes[label].xaxis.set_ticks_position("bottom")
-
-    for label in ['A', 'B', 'C']:
-        axes[label].yaxis.set_ticks_position('none')
-
-
-    """
-    Load data
-    """
-#     LOAD_ORIGINAL_DATA = True
-
-
-#     if LOAD_ORIGINAL_DATA:
-#         # use T=10500 simulation for spike raster plots
-#         label_spikes = '3afaec94d650c637ef8419611c3f80b3cb3ff539'
-#         # and T=100500 simulation for all other panels
-#         label = '99c0024eacc275d13f719afd59357f7d12f02b77'
-#         data_path = original_data_path
-#     else:
-#         from network_simulations import init_models
-#         from config import data_path
-#         models = init_models('Fig5')
-#         label_spikes = models[0].simulation.label
-#         label = models[1].simulation.label
-
-    # """
-    # Create MultiAreaModel instance to have access to data structures
-    # """
-    # 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
-    
-    # stationary firing rates
-    fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json')
-    with open(fn, 'r') as f:
-        pop_rates = json.load(f)
-
-    # time series of firing rates
-    rate_time_series = {}
-    for area in areas:
-        # fn = os.path.join(data_path, label,
-        #                   'Analysis',
-        #                   'rate_time_series_full',
-        #                   'rate_time_series_full_{}.npy'.format(area))
-        fn = os.path.join(data_path, label,
-                          'Analysis',
-                          'rate_time_series-{}.npy'.format(area))
-        rate_time_series[area] = np.load(fn)
-
-    # time series of firing rates convolved with a kernel
-    # rate_time_series_auto_kernel = {}
-    # for area in areas:
-    #     fn = os.path.join(data_path, label,
-    #                       'Analysis',
-    #                       'rate_time_series_auto_kernel',
-    #                       'rate_time_series_auto_kernel_{}.npy'.format(area))
-    #     rate_time_series_auto_kernel[area] = np.load(fn)
-
-    # local variance revised (LvR)
-    fn = os.path.join(data_path, label, 'Analysis', 'pop_LvR.json')
-    with open(fn, 'r') as f:
-        pop_LvR = json.load(f)
-
-    # correlation coefficients
-    # fn = os.path.join(data_path, label, 'Analysis', 'corrcoeff.json')
-    fn = os.path.join(data_path, label, 'Analysis', 'synchrony.json')
-    with open(fn, 'r') as f:
-        corrcoeff = json.load(f)
-
-    """
-    Plotting
-    """
-    # print("Raster plots")
-
-    # t_min = 3000.
-    # t_max = 3500.
-    t_min = t_sim - 500
-    t_max = t_sim
-
-    icolor = myred
-    ecolor = myblue
-
-    # frac_neurons = 0.03
-    frac_neurons = 1
-
-    for i, area in enumerate(areas):
-        ax = axes[labels[i]]
-
-        if area in spike_data:
-            n_pops = len(spike_data[area])
-            # Determine number of neurons that will be plotted for this area (for
-            # vertical offset)
-            offset = 0
-            n_to_plot = {}
-            for pop in M.structure[area]:
-                n_to_plot[pop] = int(M.N[area][pop] * frac_neurons)
-                offset = offset + n_to_plot[pop]
-            y_max = offset + 1
-            prev_pop = ''
-            yticks = []
-            yticklocs = []
-            for jj, pop in enumerate(M.structure[area]):
-                if pop[0:-1] != prev_pop:
-                    prev_pop = pop[0:-1]
-                    yticks.append('L' + population_labels[jj][0:-1])
-                    yticklocs.append(offset - 0.5 * n_to_plot[pop])
-                ind = np.where(np.logical_and(
-                    spike_data[area][pop][:, 1] <= t_max, spike_data[area][pop][:, 1] >= t_min))
-                pop_data = spike_data[area][pop][ind]
-                pop_neurons = np.unique(pop_data[:, 0])
-                neurons_to_ = np.arange(np.min(spike_data[area][pop][:, 0]), np.min(
-                    spike_data[area][pop][:, 0]) + n_to_plot[pop], 1)
-
-                if pop.find('E') > (-1):
-                    pcolor = ecolor
-                else:
-                    pcolor = icolor
-
-                for kk in range(n_to_plot[pop]):
-                    spike_times = pop_data[pop_data[:, 0] == neurons_to_[kk], 1]
-
-                    _ = ax.plot(spike_times, np.zeros(len(spike_times)) +
-                                offset - kk, '.', color=pcolor, markersize=1)
-                offset = offset - n_to_plot[pop]
-            y_min = offset
-            ax.set_xlim([t_min, t_max])
-            ax.set_ylim([y_min, y_max])
-            ax.set_yticklabels(yticks)
-            ax.set_yticks(yticklocs)
-            ax.set_xlabel('Time (s)', labelpad=-0.1)
-            ax.set_xticks([t_min, t_min + 250., t_max])
-            # ax.set_xticklabels([r'$3.$', r'$3.25$', r'$3.5$'])
-            l = t_min/1000
-            m = (t_min + t_max)/2000
-            r = t_max/1000
-            ax.set_xticklabels([f'{l:.2f}', f'{m:.2f}', f'{r:.2f}'])
-
-    # print("plotting Population rates")
-
-    rates = np.zeros((len(M.area_list), 8))
-    for i, area in enumerate(M.area_list):
-        for j, pop in enumerate(M.structure[area][::-1]):
-            rate = pop_rates[area][pop][0]
-            if rate == 0.0:
-                rate = 1e-5
-            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
-                rates[i][j + 2] = rate
-            else:
-                rates[i][j] = rate
-
-
-    rates = np.transpose(rates)
-    masked_rates = np.ma.masked_where(rates < 1e-4, rates)
-
-    ax = axes['D']
-    d = ax.boxplot(np.transpose(rates), vert=False,
-                   patch_artist=True, whis=1.5, showmeans=True)
-    set_boxplot_props(d)
-
-    ax.plot(np.mean(rates, axis=1), np.arange(
-        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
-    ax.set_yticklabels(population_labels[::-1], size=8)
-    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
-    ax.set_ylim((0., len(M.structure['V1']) + .5))
-
-    x_max = 220.
-    ax.set_xlim((-1., x_max))
-    ax.set_xlabel(r'Rate (spikes/s)', labelpad=-0.1)
-    ax.set_xticks([0., 50., 100.])
-
-    # print("plotting Synchrony")
-
-    syn = np.zeros((len(M.area_list), 8))
-    for i, area in enumerate(M.area_list):
-        for j, pop in enumerate(M.structure[area][::-1]):
-            value = corrcoeff[area][pop]
-            if value == 0.0:
-                value = 1e-5
-            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
-                syn[i][j + 2] = value
-            else:
-                syn[i][j] = value
-
-
-    syn = np.transpose(syn)
-    masked_syn = np.ma.masked_where(syn < 1e-4, syn)
-
-    ax = axes['E']
-    d = ax.boxplot(np.transpose(syn), vert=False,
-                   patch_artist=True, whis=1.5, showmeans=True)
-    set_boxplot_props(d)
-
-    ax.plot(np.mean(syn, axis=1), np.arange(
-        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
-
-    ax.set_yticklabels(population_labels[::-1], size=8)
-    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
-    ax.set_ylim((0., len(M.structure['V1']) + .5))
-    # ax.set_xticks(np.arange(0.0, 0.601, 0.2))
-    ax.set_xticks(np.arange(0.0, 10.0, 2))
-    ax.set_xlabel('Correlation coefficient', labelpad=-0.1)
-
-
-    # print("plotting Irregularity")
-
-    LvR = np.zeros((len(M.area_list), 8))
-    for i, area in enumerate(M.area_list):
-        for j, pop in enumerate(M.structure[area][::-1]):
-            value = pop_LvR[area][pop]
-            if value == 0.0:
-                value = 1e-5
-            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
-                LvR[i][j + 2] = value
-            else:
-                LvR[i][j] = value
-
-    LvR = np.transpose(LvR)
-    masked_LvR = np.ma.masked_where(LvR < 1e-4, LvR)
-
-    ax = axes['F']
-    d = ax.boxplot(np.transpose(LvR), vert=False,
-                   patch_artist=True, whis=1.5, showmeans=True)
-    set_boxplot_props(d)
-
-    ax.plot(np.mean(LvR, axis=1), np.arange(
-        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
-    ax.set_yticklabels(population_labels[::-1], size=8)
-    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
-    ax.set_ylim((0., len(M.structure['V1']) + .5))
-
-
-    x_max = 2.9
-    ax.set_xlim((0., x_max))
-    ax.set_xlabel('Irregularity', labelpad=-0.1)
-    ax.set_xticks([0., 1., 2.])
-
-    axes['G'].spines['right'].set_color('none')
-    axes['G'].spines['left'].set_color('none')
-    axes['G'].spines['top'].set_color('none')
-    axes['G'].spines['bottom'].set_color('none')
-    axes['G'].yaxis.set_ticks_position("none")
-    axes['G'].xaxis.set_ticks_position("none")
-    axes['G'].set_xticks([])
-    axes['G'].set_yticks([])
-
-
-    # print("Plotting rate time series")
-    pos = axes['G'].get_position()
-    ax = []
-    h = pos.y1 - pos.y0
-    w = pos.x1 - pos.x0
-    ax.append(pl.axes([pos.x0, pos.y0, w, 0.28 * h]))
-    ax.append(pl.axes([pos.x0, pos.y0 + 0.33 * h, w, 0.28 * h]))
-    ax.append(pl.axes([pos.x0, pos.y0 + 0.67 * h, w, 0.28 * h]))
-
-    colors = ['0.5', '0.3', '0.0']
-
-    # t_min = 500.
-    # t_max = 10500.
-    t_min = 500.
-    t_max = t_sim
-    # time = np.arange(500., t_max)
-    time = np.arange(500, t_max)
-    for i, area in enumerate(areas[::-1]):
-        ax[i].spines['right'].set_color('none')
-        ax[i].spines['top'].set_color('none')
-        ax[i].yaxis.set_ticks_position("left")
-        ax[i].xaxis.set_ticks_position("none")
-
-        binned_spikes = rate_time_series[area][np.where(
-            np.logical_and(time >= t_min, time < t_max))]
-        ax[i].plot(time, binned_spikes, color=colors[0], label=area)
-        # rate = rate_time_series_auto_kernel[area]
-        rate = rate_time_series[area]
-        ax[i].plot(time, rate, color=colors[2], label=area)
-        ax[i].set_xlim((t_min, t_max))
-
-        ax[i].text(0.8, 0.7, area, transform=ax[i].transAxes)
-
-        if i > 0:
-            ax[i].spines['bottom'].set_color('none')
-            ax[i].set_xticks([])
-            ax[i].set_yticks([0., 30.])
-        else:
-            # ax[i].set_xticks([1000., 5000., 10000.])
-            ax[i].set_xticks([t_min, (t_min+t_max)/2, t_max])
-            l = t_min/1000
-            m = (t_min + t_max)/2000
-            r = t_max/1000
-            # ax[i].set_xticklabels([r'$1.$', r'$5.$', r'$10.$'])
-            ax[i].set_xticklabels([f'{l:.2f}', f'{m:.2f}', f'{r:.2f}'])
-            ax[i].set_yticks([0., 5.])
-        if i == 1:
-            ax[i].set_ylabel(r'Rate (spikes/s)')
-
-    ax[0].set_xlabel('Time (s)', labelpad=-0.05)
-
-    fig.subplots_adjust(left=0.05, right=0.95, top=0.95,
-                        bottom=0.075, wspace=1., hspace=.5)
-
-    # pl.savefig('Fig5_ground_state.eps')
\ No newline at end of file
diff --git a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_instaneus_and_mean_rate-checkpoint.py b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_instaneous_and_mean_rate-checkpoint.py
similarity index 86%
rename from figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_instaneus_and_mean_rate-checkpoint.py
rename to figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_instaneous_and_mean_rate-checkpoint.py
index 1d223e89cd4911327489812e8eb16008b38eae16..065bea717950fbb888a9f7ca2d27a04932f9694b 100644
--- a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_instaneus_and_mean_rate-checkpoint.py
+++ b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_instaneous_and_mean_rate-checkpoint.py
@@ -3,13 +3,14 @@ def plot_instan_mean_firing_rate(M):
     data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + "-spikes-1-0.dat", skiprows=3)
     tsteps, spikecount = np.unique(data[:,1], return_counts=True)
     firing_rate = spikecount / M.simulation.params['dt'] * 1e3 / np.sum(M.N_vec)
-
+    
+    # visualize calculate instantaneous and mean firing rates
     ax = pl.subplot()
     ax.plot(tsteps, rate)
     ax.plot(tsteps, np.average(rate)*np.ones(len(tsteps)), label='mean')
     ax.set_title('Instantaneous and mean firing rate across all populations')
     ax.set_xlabel('time (ms)')
     ax.set_ylabel('firing rate (spikes / s)')
-    ax.set_xlim(0, sim_params['t_sim'])
+    ax.set_xlim(0, M.simulation['t_sim'])
     ax.set_ylim(0, 50)
     ax.legend()
\ No newline at end of file
diff --git a/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_resting_state-checkpoint.py b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_resting_state-checkpoint.py
new file mode 100644
index 0000000000000000000000000000000000000000..71410c0fd207c1d799c9a9afd69c826226d2e002
--- /dev/null
+++ b/figures/MAM2EBRAINS/.ipynb_checkpoints/M2E_visualize_resting_state-checkpoint.py
@@ -0,0 +1,432 @@
+import json
+import numpy as np
+import os
+
+import sys
+sys.path.append('./figures/Schmidt2018_dyn')
+
+from helpers import original_data_path, population_labels
+from multiarea_model import MultiAreaModel
+from plotcolors import myred, myblue
+
+import matplotlib.pyplot as pl
+from matplotlib import gridspec
+# from matplotlib import rc_file
+# rc_file('plotstyle.rc')
+
+icolor = myred
+ecolor = myblue
+
+# label_spikes = M.simulation.label
+label = M.simulation.label
+
+def set_boxplot_props(d):
+    for i in range(len(d['boxes'])):
+        if i % 2 == 0:
+            d['boxes'][i].set_facecolor(icolor)
+            d['boxes'][i].set_color(icolor)
+        else:
+            d['boxes'][i].set_facecolor(ecolor)
+            d['boxes'][i].set_color(ecolor)
+    pl.setp(d['whiskers'], color='k')
+    pl.setp(d['fliers'], color='k', markerfacecolor='k', marker='+')
+    pl.setp(d['medians'], color='none')
+    pl.setp(d['caps'], color='k')
+    pl.setp(d['means'], marker='x', color='k',
+            markerfacecolor='k', markeredgecolor='k', markersize=3.)
+
+def plot_resting_state(M, A, label_spikes, data_path):
+    t_sim = M.simulation.params["t_sim"]
+    
+    """
+    Figure layout
+    """
+
+    nrows = 4
+    ncols = 4
+    # width = 7.0866
+    width = 10
+    panel_wh_ratio = 0.7 * (1. + np.sqrt(5)) / 2.  # golden ratio
+
+    height = width / panel_wh_ratio * float(nrows) / ncols
+    pl.rcParams['figure.figsize'] = (width, height)
+
+
+    fig = pl.figure()
+    axes = {}
+
+    gs1 = gridspec.GridSpec(1, 3)
+    gs1.update(left=0.06, right=0.72, top=0.95, wspace=0.4, bottom=0.35)
+    # axes['A'] = pl.subplot(gs1[:-1, :1])
+    # axes['B'] = pl.subplot(gs1[:-1, 1:2])
+    # axes['C'] = pl.subplot(gs1[:-1, 2:])
+    axes['A'] = pl.subplot(gs1[:1, :1])
+    axes['B'] = pl.subplot(gs1[:1, 1:2])
+    axes['C'] = pl.subplot(gs1[:1, 2:])
+
+    gs2 = gridspec.GridSpec(3, 1)
+    gs2.update(left=0.78, right=0.95, top=0.95, bottom=0.35)
+    axes['D'] = pl.subplot(gs2[:1, :1])
+    axes['E'] = pl.subplot(gs2[1:2, :1])
+    axes['F'] = pl.subplot(gs2[2:3, :1])
+
+
+    gs3 = gridspec.GridSpec(1, 1)
+    gs3.update(left=0.1, right=0.95, top=0.3, bottom=0.075)
+    axes['G'] = pl.subplot(gs3[:1, :1])
+
+    areas = ['V1', 'V2', 'FEF']
+
+    labels = ['A', 'B', 'C']
+    for area, label in zip(areas, labels):
+        label_pos = [-0.2, 1.01]
+        # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label + ': ' + area,
+        #         fontdict={'fontsize': 10, 'weight': 'bold',
+        #                   'horizontalalignment': 'left', 'verticalalignment':
+        #                   'bottom'}, transform=axes[label].transAxes)
+        pl.text(label_pos[0], label_pos[1], label + ': ' + area,
+                 fontdict={'fontsize': 10, 'weight': 'bold', 
+                           'horizontalalignment': 'left', 'verticalalignment': 
+                           'bottom'}, transform=axes[label].transAxes)
+
+    label = 'G'
+    label_pos = [-0.1, 0.92]
+    # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
+    #         fontdict={'fontsize': 10, 'weight': 'bold',
+    #                   'horizontalalignment': 'left', 'verticalalignment':
+    #                   'bottom'}, transform=axes[label].transAxes)
+    pl.text(label_pos[0], label_pos[1], label,
+             fontdict={'fontsize': 10, 'weight': 'bold', 
+                       'horizontalalignment': 'left', 'verticalalignment': 
+                       'bottom'}, transform=axes[label].transAxes)
+
+    labels = ['E', 'D', 'F']
+    for label in labels:
+        label_pos = [-0.2, 1.05]
+        # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
+        #         fontdict={'fontsize': 10, 'weight': 'bold',
+        #                   'horizontalalignment': 'left', 'verticalalignment':
+        #                   'bottom'}, transform=axes[label].transAxes)
+        pl.text(label_pos[0], label_pos[1], label,
+             fontdict={'fontsize': 10, 'weight': 'bold', 
+                       'horizontalalignment': 'left', 'verticalalignment': 
+                       'bottom'}, transform=axes[label].transAxes)
+        
+
+    labels = ['A', 'B', 'C', 'D', 'E', 'F']
+
+    for label in labels:
+        axes[label].spines['right'].set_color('none')
+        axes[label].spines['top'].set_color('none')
+        axes[label].yaxis.set_ticks_position("left")
+        axes[label].xaxis.set_ticks_position("bottom")
+
+    for label in ['A', 'B', 'C']:
+        axes[label].yaxis.set_ticks_position('none')
+
+
+    """
+    Load data
+    """
+#     LOAD_ORIGINAL_DATA = True
+
+
+#     if LOAD_ORIGINAL_DATA:
+#         # use T=10500 simulation for spike raster plots
+#         label_spikes = '3afaec94d650c637ef8419611c3f80b3cb3ff539'
+#         # and T=100500 simulation for all other panels
+#         label = '99c0024eacc275d13f719afd59357f7d12f02b77'
+#         data_path = original_data_path
+#     else:
+#         from network_simulations import init_models
+#         from config import data_path
+#         models = init_models('Fig5')
+#         label_spikes = models[0].simulation.label
+#         label = models[1].simulation.label
+
+    # """
+    # Create MultiAreaModel instance to have access to data structures
+    # """
+    # 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
+    
+    # stationary firing rates
+    fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json')
+    with open(fn, 'r') as f:
+        pop_rates = json.load(f)
+
+    # time series of firing rates
+    rate_time_series = {}
+    for area in areas:
+        # fn = os.path.join(data_path, label,
+        #                   'Analysis',
+        #                   'rate_time_series_full',
+        #                   'rate_time_series_full_{}.npy'.format(area))
+        fn = os.path.join(data_path, label,
+                          'Analysis',
+                          'rate_time_series-{}.npy'.format(area))
+        rate_time_series[area] = np.load(fn)
+
+    # time series of firing rates convolved with a kernel
+    # rate_time_series_auto_kernel = {}
+    # for area in areas:
+    #     fn = os.path.join(data_path, label,
+    #                       'Analysis',
+    #                       'rate_time_series_auto_kernel',
+    #                       'rate_time_series_auto_kernel_{}.npy'.format(area))
+    #     rate_time_series_auto_kernel[area] = np.load(fn)
+
+    # local variance revised (LvR)
+    fn = os.path.join(data_path, label, 'Analysis', 'pop_LvR.json')
+    with open(fn, 'r') as f:
+        pop_LvR = json.load(f)
+
+    # correlation coefficients
+    # fn = os.path.join(data_path, label, 'Analysis', 'corrcoeff.json')
+    fn = os.path.join(data_path, label, 'Analysis', 'synchrony.json')
+    with open(fn, 'r') as f:
+        corrcoeff = json.load(f)
+
+    """
+    Plotting
+    """
+    # print("Raster plots")
+
+    # t_min = 3000.
+    # t_max = 3500.
+    t_min = t_sim - 500
+    t_max = t_sim
+
+    icolor = myred
+    ecolor = myblue
+
+    # frac_neurons = 0.03
+    frac_neurons = 1
+
+    for i, area in enumerate(areas):
+        ax = axes[labels[i]]
+
+        if area in spike_data:
+            n_pops = len(spike_data[area])
+            # Determine number of neurons that will be plotted for this area (for
+            # vertical offset)
+            offset = 0
+            n_to_plot = {}
+            for pop in M.structure[area]:
+                n_to_plot[pop] = int(M.N[area][pop] * frac_neurons)
+                offset = offset + n_to_plot[pop]
+            y_max = offset + 1
+            prev_pop = ''
+            yticks = []
+            yticklocs = []
+            for jj, pop in enumerate(M.structure[area]):
+                if pop[0:-1] != prev_pop:
+                    prev_pop = pop[0:-1]
+                    yticks.append('L' + population_labels[jj][0:-1])
+                    yticklocs.append(offset - 0.5 * n_to_plot[pop])
+                ind = np.where(np.logical_and(
+                    spike_data[area][pop][:, 1] <= t_max, spike_data[area][pop][:, 1] >= t_min))
+                pop_data = spike_data[area][pop][ind]
+                pop_neurons = np.unique(pop_data[:, 0])
+                neurons_to_ = np.arange(np.min(spike_data[area][pop][:, 0]), np.min(
+                    spike_data[area][pop][:, 0]) + n_to_plot[pop], 1)
+
+                if pop.find('E') > (-1):
+                    pcolor = ecolor
+                else:
+                    pcolor = icolor
+
+                for kk in range(n_to_plot[pop]):
+                    spike_times = pop_data[pop_data[:, 0] == neurons_to_[kk], 1]
+
+                    _ = ax.plot(spike_times, np.zeros(len(spike_times)) +
+                                offset - kk, '.', color=pcolor, markersize=1)
+                offset = offset - n_to_plot[pop]
+            y_min = offset
+            ax.set_xlim([t_min, t_max])
+            ax.set_ylim([y_min, y_max])
+            ax.set_yticklabels(yticks)
+            ax.set_yticks(yticklocs)
+            ax.set_xlabel('Time (s)', labelpad=-0.1)
+            ax.set_xticks([t_min, t_min + 250., t_max])
+            # ax.set_xticklabels([r'$3.$', r'$3.25$', r'$3.5$'])
+            l = t_min/1000
+            m = (t_min + t_max)/2000
+            r = t_max/1000
+            ax.set_xticklabels([f'{l:.2f}', f'{m:.2f}', f'{r:.2f}'])
+
+    # print("plotting Population rates")
+
+    rates = np.zeros((len(M.area_list), 8))
+    for i, area in enumerate(M.area_list):
+        for j, pop in enumerate(M.structure[area][::-1]):
+            rate = pop_rates[area][pop][0]
+            if rate == 0.0:
+                rate = 1e-5
+            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
+                rates[i][j + 2] = rate
+            else:
+                rates[i][j] = rate
+
+
+    rates = np.transpose(rates)
+    masked_rates = np.ma.masked_where(rates < 1e-4, rates)
+
+    ax = axes['D']
+    d = ax.boxplot(np.transpose(rates), vert=False,
+                   patch_artist=True, whis=1.5, showmeans=True)
+    set_boxplot_props(d)
+
+    ax.plot(np.mean(rates, axis=1), np.arange(
+        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
+    ax.set_yticklabels(population_labels[::-1], size=8)
+    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
+    ax.set_ylim((0., len(M.structure['V1']) + .5))
+
+    x_max = 220.
+    ax.set_xlim((-1., x_max))
+    ax.set_xlabel(r'Rate (spikes/s)', labelpad=-0.1)
+    ax.set_xticks([0., 50., 100.])
+
+    # print("plotting Synchrony")
+
+    syn = np.zeros((len(M.area_list), 8))
+    for i, area in enumerate(M.area_list):
+        for j, pop in enumerate(M.structure[area][::-1]):
+            value = corrcoeff[area][pop]
+            if value == 0.0:
+                value = 1e-5
+            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
+                syn[i][j + 2] = value
+            else:
+                syn[i][j] = value
+
+
+    syn = np.transpose(syn)
+    masked_syn = np.ma.masked_where(syn < 1e-4, syn)
+
+    ax = axes['E']
+    d = ax.boxplot(np.transpose(syn), vert=False,
+                   patch_artist=True, whis=1.5, showmeans=True)
+    set_boxplot_props(d)
+
+    ax.plot(np.mean(syn, axis=1), np.arange(
+        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
+
+    ax.set_yticklabels(population_labels[::-1], size=8)
+    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
+    ax.set_ylim((0., len(M.structure['V1']) + .5))
+    # ax.set_xticks(np.arange(0.0, 0.601, 0.2))
+    ax.set_xticks(np.arange(0.0, 10.0, 2))
+    ax.set_xlabel('Correlation coefficient', labelpad=-0.1)
+
+
+    # print("plotting Irregularity")
+
+    LvR = np.zeros((len(M.area_list), 8))
+    for i, area in enumerate(M.area_list):
+        for j, pop in enumerate(M.structure[area][::-1]):
+            value = pop_LvR[area][pop]
+            if value == 0.0:
+                value = 1e-5
+            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
+                LvR[i][j + 2] = value
+            else:
+                LvR[i][j] = value
+
+    LvR = np.transpose(LvR)
+    masked_LvR = np.ma.masked_where(LvR < 1e-4, LvR)
+
+    ax = axes['F']
+    d = ax.boxplot(np.transpose(LvR), vert=False,
+                   patch_artist=True, whis=1.5, showmeans=True)
+    set_boxplot_props(d)
+
+    ax.plot(np.mean(LvR, axis=1), np.arange(
+        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
+    ax.set_yticklabels(population_labels[::-1], size=8)
+    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
+    ax.set_ylim((0., len(M.structure['V1']) + .5))
+
+
+    x_max = 2.9
+    ax.set_xlim((0., x_max))
+    ax.set_xlabel('Irregularity', labelpad=-0.1)
+    ax.set_xticks([0., 1., 2.])
+
+    axes['G'].spines['right'].set_color('none')
+    axes['G'].spines['left'].set_color('none')
+    axes['G'].spines['top'].set_color('none')
+    axes['G'].spines['bottom'].set_color('none')
+    axes['G'].yaxis.set_ticks_position("none")
+    axes['G'].xaxis.set_ticks_position("none")
+    axes['G'].set_xticks([])
+    axes['G'].set_yticks([])
+
+
+    # print("Plotting rate time series")
+    pos = axes['G'].get_position()
+    ax = []
+    h = pos.y1 - pos.y0
+    w = pos.x1 - pos.x0
+    ax.append(pl.axes([pos.x0, pos.y0, w, 0.28 * h]))
+    ax.append(pl.axes([pos.x0, pos.y0 + 0.33 * h, w, 0.28 * h]))
+    ax.append(pl.axes([pos.x0, pos.y0 + 0.67 * h, w, 0.28 * h]))
+
+    colors = ['0.5', '0.3', '0.0']
+
+    # t_min = 500.
+    # t_max = 10500.
+    t_min = 500.
+    t_max = t_sim
+    # time = np.arange(500., t_max)
+    time = np.arange(500, t_max)
+    for i, area in enumerate(areas[::-1]):
+        ax[i].spines['right'].set_color('none')
+        ax[i].spines['top'].set_color('none')
+        ax[i].yaxis.set_ticks_position("left")
+        ax[i].xaxis.set_ticks_position("none")
+
+        binned_spikes = rate_time_series[area][np.where(
+            np.logical_and(time >= t_min, time < t_max))]
+        ax[i].plot(time, binned_spikes, color=colors[0], label=area)
+        # rate = rate_time_series_auto_kernel[area]
+        rate = rate_time_series[area]
+        ax[i].plot(time, rate, color=colors[2], label=area)
+        ax[i].set_xlim((t_min, t_max))
+
+        ax[i].text(0.8, 0.7, area, transform=ax[i].transAxes)
+
+        if i > 0:
+            ax[i].spines['bottom'].set_color('none')
+            ax[i].set_xticks([])
+            ax[i].set_yticks([0., 30.])
+        else:
+            # ax[i].set_xticks([1000., 5000., 10000.])
+            ax[i].set_xticks([t_min, (t_min+t_max)/2, t_max])
+            l = t_min/1000
+            m = (t_min + t_max)/2000
+            r = t_max/1000
+            # ax[i].set_xticklabels([r'$1.$', r'$5.$', r'$10.$'])
+            ax[i].set_xticklabels([f'{l:.2f}', f'{m:.2f}', f'{r:.2f}'])
+            ax[i].set_yticks([0., 5.])
+        if i == 1:
+            ax[i].set_ylabel(r'Rate (spikes/s)')
+
+    ax[0].set_xlabel('Time (s)', labelpad=-0.05)
+
+    fig.subplots_adjust(left=0.05, right=0.95, top=0.95,
+                        bottom=0.075, wspace=1., hspace=.5)
+
+    # pl.savefig('Fig5_ground_state.eps')
\ No newline at end of file
diff --git a/figures/MAM2EBRAINS/.ipynb_checkpoints/untitled-checkpoint.py b/figures/MAM2EBRAINS/.ipynb_checkpoints/untitled-checkpoint.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/figures/MAM2EBRAINS/M2E_VISUALIZATION.py b/figures/MAM2EBRAINS/M2E_VISUALIZATION.py
index b22bb9fd68ab69f79d77b40ec4498ef4b5ccc1b5..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391 100644
--- a/figures/MAM2EBRAINS/M2E_VISUALIZATION.py
+++ b/figures/MAM2EBRAINS/M2E_VISUALIZATION.py
@@ -1,452 +0,0 @@
-import json
-import numpy as np
-import os
-
-import sys
-sys.path.append('./figures/Schmidt2018_dyn')
-
-from helpers import original_data_path, population_labels
-from multiarea_model import MultiAreaModel
-from plotcolors import myred, myblue
-
-import matplotlib.pyplot as pl
-from matplotlib import gridspec
-# from matplotlib import rc_file
-# rc_file('plotstyle.rc')
-
-icolor = myred
-ecolor = myblue
-
-# label_spikes = M.simulation.label
-label = M.simulation.label
-
-from MAM2EBRAINS_LOAD_DATA import load_and_create_data
-A = load_and_create_data(M)
-
-def plot_instan_mean_firing_rate(M):
-    # load spike data and calculate instantaneous and mean firing rates
-    data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + "-spikes-1-0.dat", skiprows=3)
-    tsteps, spikecount = np.unique(data[:,1], return_counts=True)
-    firing_rate = spikecount / M.simulation.params['dt'] * 1e3 / np.sum(M.N_vec)
-    
-    # visualize calculate instantaneous and mean firing rates
-    ax = pl.subplot()
-    ax.plot(tsteps, rate)
-    ax.plot(tsteps, np.average(rate)*np.ones(len(tsteps)), label='mean')
-    ax.set_title('Instantaneous and mean firing rate across all populations')
-    ax.set_xlabel('time (ms)')
-    ax.set_ylabel('firing rate (spikes / s)')
-    ax.set_xlim(0, sim_params['t_sim'])
-    ax.set_ylim(0, 50)
-    ax.legend()
-
-def set_boxplot_props(d):
-    for i in range(len(d['boxes'])):
-        if i % 2 == 0:
-            d['boxes'][i].set_facecolor(icolor)
-            d['boxes'][i].set_color(icolor)
-        else:
-            d['boxes'][i].set_facecolor(ecolor)
-            d['boxes'][i].set_color(ecolor)
-    pl.setp(d['whiskers'], color='k')
-    pl.setp(d['fliers'], color='k', markerfacecolor='k', marker='+')
-    pl.setp(d['medians'], color='none')
-    pl.setp(d['caps'], color='k')
-    pl.setp(d['means'], marker='x', color='k',
-            markerfacecolor='k', markeredgecolor='k', markersize=3.)
-
-def plot_resting_state(M, A, label_spikes, data_path, sim_params):
-    t_sim = sim_params["t_sim"]
-    
-    """
-    Figure layout
-    """
-
-    nrows = 4
-    ncols = 4
-    # width = 7.0866
-    width = 10
-    panel_wh_ratio = 0.7 * (1. + np.sqrt(5)) / 2.  # golden ratio
-
-    height = width / panel_wh_ratio * float(nrows) / ncols
-    pl.rcParams['figure.figsize'] = (width, height)
-
-
-    fig = pl.figure()
-    axes = {}
-
-    gs1 = gridspec.GridSpec(1, 3)
-    gs1.update(left=0.06, right=0.72, top=0.95, wspace=0.4, bottom=0.35)
-    # axes['A'] = pl.subplot(gs1[:-1, :1])
-    # axes['B'] = pl.subplot(gs1[:-1, 1:2])
-    # axes['C'] = pl.subplot(gs1[:-1, 2:])
-    axes['A'] = pl.subplot(gs1[:1, :1])
-    axes['B'] = pl.subplot(gs1[:1, 1:2])
-    axes['C'] = pl.subplot(gs1[:1, 2:])
-
-    gs2 = gridspec.GridSpec(3, 1)
-    gs2.update(left=0.78, right=0.95, top=0.95, bottom=0.35)
-    axes['D'] = pl.subplot(gs2[:1, :1])
-    axes['E'] = pl.subplot(gs2[1:2, :1])
-    axes['F'] = pl.subplot(gs2[2:3, :1])
-
-
-    gs3 = gridspec.GridSpec(1, 1)
-    gs3.update(left=0.1, right=0.95, top=0.3, bottom=0.075)
-    axes['G'] = pl.subplot(gs3[:1, :1])
-
-    areas = ['V1', 'V2', 'FEF']
-
-    labels = ['A', 'B', 'C']
-    for area, label in zip(areas, labels):
-        label_pos = [-0.2, 1.01]
-        # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label + ': ' + area,
-        #         fontdict={'fontsize': 10, 'weight': 'bold',
-        #                   'horizontalalignment': 'left', 'verticalalignment':
-        #                   'bottom'}, transform=axes[label].transAxes)
-        pl.text(label_pos[0], label_pos[1], label + ': ' + area,
-                 fontdict={'fontsize': 10, 'weight': 'bold', 
-                           'horizontalalignment': 'left', 'verticalalignment': 
-                           'bottom'}, transform=axes[label].transAxes)
-
-    label = 'G'
-    label_pos = [-0.1, 0.92]
-    # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
-    #         fontdict={'fontsize': 10, 'weight': 'bold',
-    #                   'horizontalalignment': 'left', 'verticalalignment':
-    #                   'bottom'}, transform=axes[label].transAxes)
-    pl.text(label_pos[0], label_pos[1], label,
-             fontdict={'fontsize': 10, 'weight': 'bold', 
-                       'horizontalalignment': 'left', 'verticalalignment': 
-                       'bottom'}, transform=axes[label].transAxes)
-
-    labels = ['E', 'D', 'F']
-    for label in labels:
-        label_pos = [-0.2, 1.05]
-        # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
-        #         fontdict={'fontsize': 10, 'weight': 'bold',
-        #                   'horizontalalignment': 'left', 'verticalalignment':
-        #                   'bottom'}, transform=axes[label].transAxes)
-        pl.text(label_pos[0], label_pos[1], label,
-             fontdict={'fontsize': 10, 'weight': 'bold', 
-                       'horizontalalignment': 'left', 'verticalalignment': 
-                       'bottom'}, transform=axes[label].transAxes)
-        
-
-    labels = ['A', 'B', 'C', 'D', 'E', 'F']
-
-    for label in labels:
-        axes[label].spines['right'].set_color('none')
-        axes[label].spines['top'].set_color('none')
-        axes[label].yaxis.set_ticks_position("left")
-        axes[label].xaxis.set_ticks_position("bottom")
-
-    for label in ['A', 'B', 'C']:
-        axes[label].yaxis.set_ticks_position('none')
-
-
-    """
-    Load data
-    """
-#     LOAD_ORIGINAL_DATA = True
-
-
-#     if LOAD_ORIGINAL_DATA:
-#         # use T=10500 simulation for spike raster plots
-#         label_spikes = '3afaec94d650c637ef8419611c3f80b3cb3ff539'
-#         # and T=100500 simulation for all other panels
-#         label = '99c0024eacc275d13f719afd59357f7d12f02b77'
-#         data_path = original_data_path
-#     else:
-#         from network_simulations import init_models
-#         from config import data_path
-#         models = init_models('Fig5')
-#         label_spikes = models[0].simulation.label
-#         label = models[1].simulation.label
-
-    # """
-    # Create MultiAreaModel instance to have access to data structures
-    # """
-    # 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
-    
-    # stationary firing rates
-    fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json')
-    with open(fn, 'r') as f:
-        pop_rates = json.load(f)
-
-    # time series of firing rates
-    rate_time_series = {}
-    for area in areas:
-        # fn = os.path.join(data_path, label,
-        #                   'Analysis',
-        #                   'rate_time_series_full',
-        #                   'rate_time_series_full_{}.npy'.format(area))
-        fn = os.path.join(data_path, label,
-                          'Analysis',
-                          'rate_time_series-{}.npy'.format(area))
-        rate_time_series[area] = np.load(fn)
-
-    # time series of firing rates convolved with a kernel
-    # rate_time_series_auto_kernel = {}
-    # for area in areas:
-    #     fn = os.path.join(data_path, label,
-    #                       'Analysis',
-    #                       'rate_time_series_auto_kernel',
-    #                       'rate_time_series_auto_kernel_{}.npy'.format(area))
-    #     rate_time_series_auto_kernel[area] = np.load(fn)
-
-    # local variance revised (LvR)
-    fn = os.path.join(data_path, label, 'Analysis', 'pop_LvR.json')
-    with open(fn, 'r') as f:
-        pop_LvR = json.load(f)
-
-    # correlation coefficients
-    # fn = os.path.join(data_path, label, 'Analysis', 'corrcoeff.json')
-    fn = os.path.join(data_path, label, 'Analysis', 'synchrony.json')
-    with open(fn, 'r') as f:
-        corrcoeff = json.load(f)
-
-    """
-    Plotting
-    """
-    # print("Raster plots")
-
-    # t_min = 3000.
-    # t_max = 3500.
-    t_min = t_sim - 500
-    t_max = t_sim
-
-    icolor = myred
-    ecolor = myblue
-
-    # frac_neurons = 0.03
-    frac_neurons = 1
-
-    for i, area in enumerate(areas):
-        ax = axes[labels[i]]
-
-        if area in spike_data:
-            n_pops = len(spike_data[area])
-            # Determine number of neurons that will be plotted for this area (for
-            # vertical offset)
-            offset = 0
-            n_to_plot = {}
-            for pop in M.structure[area]:
-                n_to_plot[pop] = int(M.N[area][pop] * frac_neurons)
-                offset = offset + n_to_plot[pop]
-            y_max = offset + 1
-            prev_pop = ''
-            yticks = []
-            yticklocs = []
-            for jj, pop in enumerate(M.structure[area]):
-                if pop[0:-1] != prev_pop:
-                    prev_pop = pop[0:-1]
-                    yticks.append('L' + population_labels[jj][0:-1])
-                    yticklocs.append(offset - 0.5 * n_to_plot[pop])
-                ind = np.where(np.logical_and(
-                    spike_data[area][pop][:, 1] <= t_max, spike_data[area][pop][:, 1] >= t_min))
-                pop_data = spike_data[area][pop][ind]
-                pop_neurons = np.unique(pop_data[:, 0])
-                neurons_to_ = np.arange(np.min(spike_data[area][pop][:, 0]), np.min(
-                    spike_data[area][pop][:, 0]) + n_to_plot[pop], 1)
-
-                if pop.find('E') > (-1):
-                    pcolor = ecolor
-                else:
-                    pcolor = icolor
-
-                for kk in range(n_to_plot[pop]):
-                    spike_times = pop_data[pop_data[:, 0] == neurons_to_[kk], 1]
-
-                    _ = ax.plot(spike_times, np.zeros(len(spike_times)) +
-                                offset - kk, '.', color=pcolor, markersize=1)
-                offset = offset - n_to_plot[pop]
-            y_min = offset
-            ax.set_xlim([t_min, t_max])
-            ax.set_ylim([y_min, y_max])
-            ax.set_yticklabels(yticks)
-            ax.set_yticks(yticklocs)
-            ax.set_xlabel('Time (s)', labelpad=-0.1)
-            ax.set_xticks([t_min, t_min + 250., t_max])
-            # ax.set_xticklabels([r'$3.$', r'$3.25$', r'$3.5$'])
-            l = t_min/1000
-            m = (t_min + t_max)/2000
-            r = t_max/1000
-            ax.set_xticklabels([f'{l:.2f}', f'{m:.2f}', f'{r:.2f}'])
-
-    # print("plotting Population rates")
-
-    rates = np.zeros((len(M.area_list), 8))
-    for i, area in enumerate(M.area_list):
-        for j, pop in enumerate(M.structure[area][::-1]):
-            rate = pop_rates[area][pop][0]
-            if rate == 0.0:
-                rate = 1e-5
-            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
-                rates[i][j + 2] = rate
-            else:
-                rates[i][j] = rate
-
-
-    rates = np.transpose(rates)
-    masked_rates = np.ma.masked_where(rates < 1e-4, rates)
-
-    ax = axes['D']
-    d = ax.boxplot(np.transpose(rates), vert=False,
-                   patch_artist=True, whis=1.5, showmeans=True)
-    set_boxplot_props(d)
-
-    ax.plot(np.mean(rates, axis=1), np.arange(
-        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
-    ax.set_yticklabels(population_labels[::-1], size=8)
-    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
-    ax.set_ylim((0., len(M.structure['V1']) + .5))
-
-    x_max = 220.
-    ax.set_xlim((-1., x_max))
-    ax.set_xlabel(r'Rate (spikes/s)', labelpad=-0.1)
-    ax.set_xticks([0., 50., 100.])
-
-    # print("plotting Synchrony")
-
-    syn = np.zeros((len(M.area_list), 8))
-    for i, area in enumerate(M.area_list):
-        for j, pop in enumerate(M.structure[area][::-1]):
-            value = corrcoeff[area][pop]
-            if value == 0.0:
-                value = 1e-5
-            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
-                syn[i][j + 2] = value
-            else:
-                syn[i][j] = value
-
-
-    syn = np.transpose(syn)
-    masked_syn = np.ma.masked_where(syn < 1e-4, syn)
-
-    ax = axes['E']
-    d = ax.boxplot(np.transpose(syn), vert=False,
-                   patch_artist=True, whis=1.5, showmeans=True)
-    set_boxplot_props(d)
-
-    ax.plot(np.mean(syn, axis=1), np.arange(
-        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
-
-    ax.set_yticklabels(population_labels[::-1], size=8)
-    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
-    ax.set_ylim((0., len(M.structure['V1']) + .5))
-    # ax.set_xticks(np.arange(0.0, 0.601, 0.2))
-    ax.set_xticks(np.arange(0.0, 10.0, 2))
-    ax.set_xlabel('Correlation coefficient', labelpad=-0.1)
-
-
-    # print("plotting Irregularity")
-
-    LvR = np.zeros((len(M.area_list), 8))
-    for i, area in enumerate(M.area_list):
-        for j, pop in enumerate(M.structure[area][::-1]):
-            value = pop_LvR[area][pop]
-            if value == 0.0:
-                value = 1e-5
-            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
-                LvR[i][j + 2] = value
-            else:
-                LvR[i][j] = value
-
-    LvR = np.transpose(LvR)
-    masked_LvR = np.ma.masked_where(LvR < 1e-4, LvR)
-
-    ax = axes['F']
-    d = ax.boxplot(np.transpose(LvR), vert=False,
-                   patch_artist=True, whis=1.5, showmeans=True)
-    set_boxplot_props(d)
-
-    ax.plot(np.mean(LvR, axis=1), np.arange(
-        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
-    ax.set_yticklabels(population_labels[::-1], size=8)
-    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
-    ax.set_ylim((0., len(M.structure['V1']) + .5))
-
-
-    x_max = 2.9
-    ax.set_xlim((0., x_max))
-    ax.set_xlabel('Irregularity', labelpad=-0.1)
-    ax.set_xticks([0., 1., 2.])
-
-    axes['G'].spines['right'].set_color('none')
-    axes['G'].spines['left'].set_color('none')
-    axes['G'].spines['top'].set_color('none')
-    axes['G'].spines['bottom'].set_color('none')
-    axes['G'].yaxis.set_ticks_position("none")
-    axes['G'].xaxis.set_ticks_position("none")
-    axes['G'].set_xticks([])
-    axes['G'].set_yticks([])
-
-
-    # print("Plotting rate time series")
-    pos = axes['G'].get_position()
-    ax = []
-    h = pos.y1 - pos.y0
-    w = pos.x1 - pos.x0
-    ax.append(pl.axes([pos.x0, pos.y0, w, 0.28 * h]))
-    ax.append(pl.axes([pos.x0, pos.y0 + 0.33 * h, w, 0.28 * h]))
-    ax.append(pl.axes([pos.x0, pos.y0 + 0.67 * h, w, 0.28 * h]))
-
-    colors = ['0.5', '0.3', '0.0']
-
-    # t_min = 500.
-    # t_max = 10500.
-    t_min = 500.
-    t_max = t_sim
-    # time = np.arange(500., t_max)
-    time = np.arange(500, t_max)
-    for i, area in enumerate(areas[::-1]):
-        ax[i].spines['right'].set_color('none')
-        ax[i].spines['top'].set_color('none')
-        ax[i].yaxis.set_ticks_position("left")
-        ax[i].xaxis.set_ticks_position("none")
-
-        binned_spikes = rate_time_series[area][np.where(
-            np.logical_and(time >= t_min, time < t_max))]
-        ax[i].plot(time, binned_spikes, color=colors[0], label=area)
-        # rate = rate_time_series_auto_kernel[area]
-        rate = rate_time_series[area]
-        ax[i].plot(time, rate, color=colors[2], label=area)
-        ax[i].set_xlim((t_min, t_max))
-
-        ax[i].text(0.8, 0.7, area, transform=ax[i].transAxes)
-
-        if i > 0:
-            ax[i].spines['bottom'].set_color('none')
-            ax[i].set_xticks([])
-            ax[i].set_yticks([0., 30.])
-        else:
-            # ax[i].set_xticks([1000., 5000., 10000.])
-            ax[i].set_xticks([t_min, (t_min+t_max)/2, t_max])
-            l = t_min/1000
-            m = (t_min + t_max)/2000
-            r = t_max/1000
-            # ax[i].set_xticklabels([r'$1.$', r'$5.$', r'$10.$'])
-            ax[i].set_xticklabels([f'{l:.2f}', f'{m:.2f}', f'{r:.2f}'])
-            ax[i].set_yticks([0., 5.])
-        if i == 1:
-            ax[i].set_ylabel(r'Rate (spikes/s)')
-
-    ax[0].set_xlabel('Time (s)', labelpad=-0.05)
-
-    fig.subplots_adjust(left=0.05, right=0.95, top=0.95,
-                        bottom=0.075, wspace=1., hspace=.5)
-
-    # pl.savefig('Fig5_ground_state.eps')
\ No newline at end of file
diff --git a/figures/MAM2EBRAINS/M2E_visualize_instaneous_and_mean_rate.py b/figures/MAM2EBRAINS/M2E_visualize_instaneous_and_mean_rate.py
new file mode 100644
index 0000000000000000000000000000000000000000..065bea717950fbb888a9f7ca2d27a04932f9694b
--- /dev/null
+++ b/figures/MAM2EBRAINS/M2E_visualize_instaneous_and_mean_rate.py
@@ -0,0 +1,16 @@
+def plot_instan_mean_firing_rate(M):
+    # load spike data and calculate instantaneous and mean firing rates
+    data = np.loadtxt(M.simulation.data_dir + '/recordings/' + M.simulation.label + "-spikes-1-0.dat", skiprows=3)
+    tsteps, spikecount = np.unique(data[:,1], return_counts=True)
+    firing_rate = spikecount / M.simulation.params['dt'] * 1e3 / np.sum(M.N_vec)
+    
+    # visualize calculate instantaneous and mean firing rates
+    ax = pl.subplot()
+    ax.plot(tsteps, rate)
+    ax.plot(tsteps, np.average(rate)*np.ones(len(tsteps)), label='mean')
+    ax.set_title('Instantaneous and mean firing rate across all populations')
+    ax.set_xlabel('time (ms)')
+    ax.set_ylabel('firing rate (spikes / s)')
+    ax.set_xlim(0, M.simulation['t_sim'])
+    ax.set_ylim(0, 50)
+    ax.legend()
\ No newline at end of file
diff --git a/figures/MAM2EBRAINS/M2E_visualize_instaneus_and_mean_rate.py b/figures/MAM2EBRAINS/M2E_visualize_instaneus_and_mean_rate.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/figures/MAM2EBRAINS/M2E_visualize_resting_state.py b/figures/MAM2EBRAINS/M2E_visualize_resting_state.py
new file mode 100644
index 0000000000000000000000000000000000000000..71410c0fd207c1d799c9a9afd69c826226d2e002
--- /dev/null
+++ b/figures/MAM2EBRAINS/M2E_visualize_resting_state.py
@@ -0,0 +1,432 @@
+import json
+import numpy as np
+import os
+
+import sys
+sys.path.append('./figures/Schmidt2018_dyn')
+
+from helpers import original_data_path, population_labels
+from multiarea_model import MultiAreaModel
+from plotcolors import myred, myblue
+
+import matplotlib.pyplot as pl
+from matplotlib import gridspec
+# from matplotlib import rc_file
+# rc_file('plotstyle.rc')
+
+icolor = myred
+ecolor = myblue
+
+# label_spikes = M.simulation.label
+label = M.simulation.label
+
+def set_boxplot_props(d):
+    for i in range(len(d['boxes'])):
+        if i % 2 == 0:
+            d['boxes'][i].set_facecolor(icolor)
+            d['boxes'][i].set_color(icolor)
+        else:
+            d['boxes'][i].set_facecolor(ecolor)
+            d['boxes'][i].set_color(ecolor)
+    pl.setp(d['whiskers'], color='k')
+    pl.setp(d['fliers'], color='k', markerfacecolor='k', marker='+')
+    pl.setp(d['medians'], color='none')
+    pl.setp(d['caps'], color='k')
+    pl.setp(d['means'], marker='x', color='k',
+            markerfacecolor='k', markeredgecolor='k', markersize=3.)
+
+def plot_resting_state(M, A, label_spikes, data_path):
+    t_sim = M.simulation.params["t_sim"]
+    
+    """
+    Figure layout
+    """
+
+    nrows = 4
+    ncols = 4
+    # width = 7.0866
+    width = 10
+    panel_wh_ratio = 0.7 * (1. + np.sqrt(5)) / 2.  # golden ratio
+
+    height = width / panel_wh_ratio * float(nrows) / ncols
+    pl.rcParams['figure.figsize'] = (width, height)
+
+
+    fig = pl.figure()
+    axes = {}
+
+    gs1 = gridspec.GridSpec(1, 3)
+    gs1.update(left=0.06, right=0.72, top=0.95, wspace=0.4, bottom=0.35)
+    # axes['A'] = pl.subplot(gs1[:-1, :1])
+    # axes['B'] = pl.subplot(gs1[:-1, 1:2])
+    # axes['C'] = pl.subplot(gs1[:-1, 2:])
+    axes['A'] = pl.subplot(gs1[:1, :1])
+    axes['B'] = pl.subplot(gs1[:1, 1:2])
+    axes['C'] = pl.subplot(gs1[:1, 2:])
+
+    gs2 = gridspec.GridSpec(3, 1)
+    gs2.update(left=0.78, right=0.95, top=0.95, bottom=0.35)
+    axes['D'] = pl.subplot(gs2[:1, :1])
+    axes['E'] = pl.subplot(gs2[1:2, :1])
+    axes['F'] = pl.subplot(gs2[2:3, :1])
+
+
+    gs3 = gridspec.GridSpec(1, 1)
+    gs3.update(left=0.1, right=0.95, top=0.3, bottom=0.075)
+    axes['G'] = pl.subplot(gs3[:1, :1])
+
+    areas = ['V1', 'V2', 'FEF']
+
+    labels = ['A', 'B', 'C']
+    for area, label in zip(areas, labels):
+        label_pos = [-0.2, 1.01]
+        # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label + ': ' + area,
+        #         fontdict={'fontsize': 10, 'weight': 'bold',
+        #                   'horizontalalignment': 'left', 'verticalalignment':
+        #                   'bottom'}, transform=axes[label].transAxes)
+        pl.text(label_pos[0], label_pos[1], label + ': ' + area,
+                 fontdict={'fontsize': 10, 'weight': 'bold', 
+                           'horizontalalignment': 'left', 'verticalalignment': 
+                           'bottom'}, transform=axes[label].transAxes)
+
+    label = 'G'
+    label_pos = [-0.1, 0.92]
+    # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
+    #         fontdict={'fontsize': 10, 'weight': 'bold',
+    #                   'horizontalalignment': 'left', 'verticalalignment':
+    #                   'bottom'}, transform=axes[label].transAxes)
+    pl.text(label_pos[0], label_pos[1], label,
+             fontdict={'fontsize': 10, 'weight': 'bold', 
+                       'horizontalalignment': 'left', 'verticalalignment': 
+                       'bottom'}, transform=axes[label].transAxes)
+
+    labels = ['E', 'D', 'F']
+    for label in labels:
+        label_pos = [-0.2, 1.05]
+        # pl.text(label_pos[0], label_pos[1], r'\bfseries{}' + label,
+        #         fontdict={'fontsize': 10, 'weight': 'bold',
+        #                   'horizontalalignment': 'left', 'verticalalignment':
+        #                   'bottom'}, transform=axes[label].transAxes)
+        pl.text(label_pos[0], label_pos[1], label,
+             fontdict={'fontsize': 10, 'weight': 'bold', 
+                       'horizontalalignment': 'left', 'verticalalignment': 
+                       'bottom'}, transform=axes[label].transAxes)
+        
+
+    labels = ['A', 'B', 'C', 'D', 'E', 'F']
+
+    for label in labels:
+        axes[label].spines['right'].set_color('none')
+        axes[label].spines['top'].set_color('none')
+        axes[label].yaxis.set_ticks_position("left")
+        axes[label].xaxis.set_ticks_position("bottom")
+
+    for label in ['A', 'B', 'C']:
+        axes[label].yaxis.set_ticks_position('none')
+
+
+    """
+    Load data
+    """
+#     LOAD_ORIGINAL_DATA = True
+
+
+#     if LOAD_ORIGINAL_DATA:
+#         # use T=10500 simulation for spike raster plots
+#         label_spikes = '3afaec94d650c637ef8419611c3f80b3cb3ff539'
+#         # and T=100500 simulation for all other panels
+#         label = '99c0024eacc275d13f719afd59357f7d12f02b77'
+#         data_path = original_data_path
+#     else:
+#         from network_simulations import init_models
+#         from config import data_path
+#         models = init_models('Fig5')
+#         label_spikes = models[0].simulation.label
+#         label = models[1].simulation.label
+
+    # """
+    # Create MultiAreaModel instance to have access to data structures
+    # """
+    # 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
+    
+    # stationary firing rates
+    fn = os.path.join(data_path, label, 'Analysis', 'pop_rates.json')
+    with open(fn, 'r') as f:
+        pop_rates = json.load(f)
+
+    # time series of firing rates
+    rate_time_series = {}
+    for area in areas:
+        # fn = os.path.join(data_path, label,
+        #                   'Analysis',
+        #                   'rate_time_series_full',
+        #                   'rate_time_series_full_{}.npy'.format(area))
+        fn = os.path.join(data_path, label,
+                          'Analysis',
+                          'rate_time_series-{}.npy'.format(area))
+        rate_time_series[area] = np.load(fn)
+
+    # time series of firing rates convolved with a kernel
+    # rate_time_series_auto_kernel = {}
+    # for area in areas:
+    #     fn = os.path.join(data_path, label,
+    #                       'Analysis',
+    #                       'rate_time_series_auto_kernel',
+    #                       'rate_time_series_auto_kernel_{}.npy'.format(area))
+    #     rate_time_series_auto_kernel[area] = np.load(fn)
+
+    # local variance revised (LvR)
+    fn = os.path.join(data_path, label, 'Analysis', 'pop_LvR.json')
+    with open(fn, 'r') as f:
+        pop_LvR = json.load(f)
+
+    # correlation coefficients
+    # fn = os.path.join(data_path, label, 'Analysis', 'corrcoeff.json')
+    fn = os.path.join(data_path, label, 'Analysis', 'synchrony.json')
+    with open(fn, 'r') as f:
+        corrcoeff = json.load(f)
+
+    """
+    Plotting
+    """
+    # print("Raster plots")
+
+    # t_min = 3000.
+    # t_max = 3500.
+    t_min = t_sim - 500
+    t_max = t_sim
+
+    icolor = myred
+    ecolor = myblue
+
+    # frac_neurons = 0.03
+    frac_neurons = 1
+
+    for i, area in enumerate(areas):
+        ax = axes[labels[i]]
+
+        if area in spike_data:
+            n_pops = len(spike_data[area])
+            # Determine number of neurons that will be plotted for this area (for
+            # vertical offset)
+            offset = 0
+            n_to_plot = {}
+            for pop in M.structure[area]:
+                n_to_plot[pop] = int(M.N[area][pop] * frac_neurons)
+                offset = offset + n_to_plot[pop]
+            y_max = offset + 1
+            prev_pop = ''
+            yticks = []
+            yticklocs = []
+            for jj, pop in enumerate(M.structure[area]):
+                if pop[0:-1] != prev_pop:
+                    prev_pop = pop[0:-1]
+                    yticks.append('L' + population_labels[jj][0:-1])
+                    yticklocs.append(offset - 0.5 * n_to_plot[pop])
+                ind = np.where(np.logical_and(
+                    spike_data[area][pop][:, 1] <= t_max, spike_data[area][pop][:, 1] >= t_min))
+                pop_data = spike_data[area][pop][ind]
+                pop_neurons = np.unique(pop_data[:, 0])
+                neurons_to_ = np.arange(np.min(spike_data[area][pop][:, 0]), np.min(
+                    spike_data[area][pop][:, 0]) + n_to_plot[pop], 1)
+
+                if pop.find('E') > (-1):
+                    pcolor = ecolor
+                else:
+                    pcolor = icolor
+
+                for kk in range(n_to_plot[pop]):
+                    spike_times = pop_data[pop_data[:, 0] == neurons_to_[kk], 1]
+
+                    _ = ax.plot(spike_times, np.zeros(len(spike_times)) +
+                                offset - kk, '.', color=pcolor, markersize=1)
+                offset = offset - n_to_plot[pop]
+            y_min = offset
+            ax.set_xlim([t_min, t_max])
+            ax.set_ylim([y_min, y_max])
+            ax.set_yticklabels(yticks)
+            ax.set_yticks(yticklocs)
+            ax.set_xlabel('Time (s)', labelpad=-0.1)
+            ax.set_xticks([t_min, t_min + 250., t_max])
+            # ax.set_xticklabels([r'$3.$', r'$3.25$', r'$3.5$'])
+            l = t_min/1000
+            m = (t_min + t_max)/2000
+            r = t_max/1000
+            ax.set_xticklabels([f'{l:.2f}', f'{m:.2f}', f'{r:.2f}'])
+
+    # print("plotting Population rates")
+
+    rates = np.zeros((len(M.area_list), 8))
+    for i, area in enumerate(M.area_list):
+        for j, pop in enumerate(M.structure[area][::-1]):
+            rate = pop_rates[area][pop][0]
+            if rate == 0.0:
+                rate = 1e-5
+            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
+                rates[i][j + 2] = rate
+            else:
+                rates[i][j] = rate
+
+
+    rates = np.transpose(rates)
+    masked_rates = np.ma.masked_where(rates < 1e-4, rates)
+
+    ax = axes['D']
+    d = ax.boxplot(np.transpose(rates), vert=False,
+                   patch_artist=True, whis=1.5, showmeans=True)
+    set_boxplot_props(d)
+
+    ax.plot(np.mean(rates, axis=1), np.arange(
+        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
+    ax.set_yticklabels(population_labels[::-1], size=8)
+    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
+    ax.set_ylim((0., len(M.structure['V1']) + .5))
+
+    x_max = 220.
+    ax.set_xlim((-1., x_max))
+    ax.set_xlabel(r'Rate (spikes/s)', labelpad=-0.1)
+    ax.set_xticks([0., 50., 100.])
+
+    # print("plotting Synchrony")
+
+    syn = np.zeros((len(M.area_list), 8))
+    for i, area in enumerate(M.area_list):
+        for j, pop in enumerate(M.structure[area][::-1]):
+            value = corrcoeff[area][pop]
+            if value == 0.0:
+                value = 1e-5
+            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
+                syn[i][j + 2] = value
+            else:
+                syn[i][j] = value
+
+
+    syn = np.transpose(syn)
+    masked_syn = np.ma.masked_where(syn < 1e-4, syn)
+
+    ax = axes['E']
+    d = ax.boxplot(np.transpose(syn), vert=False,
+                   patch_artist=True, whis=1.5, showmeans=True)
+    set_boxplot_props(d)
+
+    ax.plot(np.mean(syn, axis=1), np.arange(
+        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
+
+    ax.set_yticklabels(population_labels[::-1], size=8)
+    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
+    ax.set_ylim((0., len(M.structure['V1']) + .5))
+    # ax.set_xticks(np.arange(0.0, 0.601, 0.2))
+    ax.set_xticks(np.arange(0.0, 10.0, 2))
+    ax.set_xlabel('Correlation coefficient', labelpad=-0.1)
+
+
+    # print("plotting Irregularity")
+
+    LvR = np.zeros((len(M.area_list), 8))
+    for i, area in enumerate(M.area_list):
+        for j, pop in enumerate(M.structure[area][::-1]):
+            value = pop_LvR[area][pop]
+            if value == 0.0:
+                value = 1e-5
+            if area == 'TH' and j > 3:  # To account for missing layer 4 in TH
+                LvR[i][j + 2] = value
+            else:
+                LvR[i][j] = value
+
+    LvR = np.transpose(LvR)
+    masked_LvR = np.ma.masked_where(LvR < 1e-4, LvR)
+
+    ax = axes['F']
+    d = ax.boxplot(np.transpose(LvR), vert=False,
+                   patch_artist=True, whis=1.5, showmeans=True)
+    set_boxplot_props(d)
+
+    ax.plot(np.mean(LvR, axis=1), np.arange(
+        1., len(M.structure['V1']) + 1., 1.), 'x', color='k', markersize=3)
+    ax.set_yticklabels(population_labels[::-1], size=8)
+    ax.set_yticks(np.arange(1., len(M.structure['V1']) + 1., 1.))
+    ax.set_ylim((0., len(M.structure['V1']) + .5))
+
+
+    x_max = 2.9
+    ax.set_xlim((0., x_max))
+    ax.set_xlabel('Irregularity', labelpad=-0.1)
+    ax.set_xticks([0., 1., 2.])
+
+    axes['G'].spines['right'].set_color('none')
+    axes['G'].spines['left'].set_color('none')
+    axes['G'].spines['top'].set_color('none')
+    axes['G'].spines['bottom'].set_color('none')
+    axes['G'].yaxis.set_ticks_position("none")
+    axes['G'].xaxis.set_ticks_position("none")
+    axes['G'].set_xticks([])
+    axes['G'].set_yticks([])
+
+
+    # print("Plotting rate time series")
+    pos = axes['G'].get_position()
+    ax = []
+    h = pos.y1 - pos.y0
+    w = pos.x1 - pos.x0
+    ax.append(pl.axes([pos.x0, pos.y0, w, 0.28 * h]))
+    ax.append(pl.axes([pos.x0, pos.y0 + 0.33 * h, w, 0.28 * h]))
+    ax.append(pl.axes([pos.x0, pos.y0 + 0.67 * h, w, 0.28 * h]))
+
+    colors = ['0.5', '0.3', '0.0']
+
+    # t_min = 500.
+    # t_max = 10500.
+    t_min = 500.
+    t_max = t_sim
+    # time = np.arange(500., t_max)
+    time = np.arange(500, t_max)
+    for i, area in enumerate(areas[::-1]):
+        ax[i].spines['right'].set_color('none')
+        ax[i].spines['top'].set_color('none')
+        ax[i].yaxis.set_ticks_position("left")
+        ax[i].xaxis.set_ticks_position("none")
+
+        binned_spikes = rate_time_series[area][np.where(
+            np.logical_and(time >= t_min, time < t_max))]
+        ax[i].plot(time, binned_spikes, color=colors[0], label=area)
+        # rate = rate_time_series_auto_kernel[area]
+        rate = rate_time_series[area]
+        ax[i].plot(time, rate, color=colors[2], label=area)
+        ax[i].set_xlim((t_min, t_max))
+
+        ax[i].text(0.8, 0.7, area, transform=ax[i].transAxes)
+
+        if i > 0:
+            ax[i].spines['bottom'].set_color('none')
+            ax[i].set_xticks([])
+            ax[i].set_yticks([0., 30.])
+        else:
+            # ax[i].set_xticks([1000., 5000., 10000.])
+            ax[i].set_xticks([t_min, (t_min+t_max)/2, t_max])
+            l = t_min/1000
+            m = (t_min + t_max)/2000
+            r = t_max/1000
+            # ax[i].set_xticklabels([r'$1.$', r'$5.$', r'$10.$'])
+            ax[i].set_xticklabels([f'{l:.2f}', f'{m:.2f}', f'{r:.2f}'])
+            ax[i].set_yticks([0., 5.])
+        if i == 1:
+            ax[i].set_ylabel(r'Rate (spikes/s)')
+
+    ax[0].set_xlabel('Time (s)', labelpad=-0.05)
+
+    fig.subplots_adjust(left=0.05, right=0.95, top=0.95,
+                        bottom=0.075, wspace=1., hspace=.5)
+
+    # pl.savefig('Fig5_ground_state.eps')
\ No newline at end of file
diff --git a/figures/MAM2EBRAINS/untitled.py b/figures/MAM2EBRAINS/untitled.py
deleted file mode 100644
index e69de29bb2d1d6434b8b29ae775ad8c2e48c5391..0000000000000000000000000000000000000000
diff --git a/multi-area-model.ipynb b/multi-area-model.ipynb
index baef6de81efd308c240de5fdd9eb76f20ff46d8b..e766a472f4359c459c13afe732002d489e9556b0 100644
--- a/multi-area-model.ipynb
+++ b/multi-area-model.ipynb
@@ -27,11 +27,10 @@
     "    * [2.2. Predict firing rates from theory](#section_2_2)\n",
     "    * [2.3. Extract and visualize interareal connectivity](#section_2_3)\n",
     "    * [2.4. Run a simulation](#section_2_4)\n",
-    "* [S3. Data Loading and Processing](#section_3)\n",
-    "* [S4. Simulation Results Visualization](#section_4) \n",
-    "    * [4.1. Instantaneous and mean firing rate across all populations](#section_4_1)\n",
-    "    * [4.2 Resting state plots](#section_4_2)\n",
-    "    * [4.3 Time-averaged population rates](#section_4_3)"
+    "* [S3. Simulation Results Visualization](#section_3) \n",
+    "    * [3.1. Instantaneous and mean firing rate across all populations](#section_3_1)\n",
+    "    * [3.2. Resting state plots](#section_3_2)\n",
+    "    * [3.3. Time-averaged population rates](#section_3_3)"
    ]
   },
   {
@@ -46,7 +45,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 37,
+   "execution_count": 49,
    "id": "9d6cc7d9-3110-4d96-9f9a-9ec7dee6d145",
    "metadata": {},
    "outputs": [],
@@ -64,7 +63,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 38,
+   "execution_count": 53,
    "id": "96517739",
    "metadata": {
     "tags": []
@@ -77,25 +76,17 @@
     "import nest\n",
     "import json\n",
     "import sys\n",
-    "from io import StringIO\n",
+    "from IPython.display import display, HTML\n",
     "\n",
     "from multiarea_model import MultiAreaModel\n",
     "from config import base_path, data_path\n",
     "\n",
-    "sys.path.append('./figures/MAM2EBRAINS')\n",
-    "\n",
-    "# # Ignore and don't display warning messages\n",
-    "# import warnings\n",
-    "# warnings.filterwarnings(\"ignore\")\n",
-    "\n",
-    "# # Redirect stdout to a dummy stream\n",
-    "# original_stdout = sys.stdout\n",
-    "# sys.stdout = StringIO()"
+    "sys.path.append('./figures/MAM2EBRAINS')"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 39,
+   "execution_count": 54,
    "id": "7e07b0d0",
    "metadata": {
     "tags": []
@@ -108,7 +99,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 40,
+   "execution_count": 55,
    "id": "1d440c07-9b69-4e52-8573-26b13493bc5a",
    "metadata": {
     "tags": []
@@ -132,7 +123,6 @@
    ],
    "source": [
     "# Jupyter notebook display format setting\n",
-    "from IPython.display import display, HTML\n",
     "style = \"\"\"\n",
     "<style>\n",
     "table {float:left}\n",
@@ -199,7 +189,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 41,
+   "execution_count": 56,
    "id": "60265d52",
    "metadata": {},
    "outputs": [],
@@ -236,7 +226,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 42,
+   "execution_count": 57,
    "id": "6e4bed8d",
    "metadata": {},
    "outputs": [],
@@ -319,7 +309,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 43,
+   "execution_count": 58,
    "id": "ab25f9f8",
    "metadata": {},
    "outputs": [
@@ -351,7 +341,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 44,
+   "execution_count": 59,
    "id": "6a7ddf0e",
    "metadata": {},
    "outputs": [],
@@ -379,7 +369,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 45,
+   "execution_count": 60,
    "id": "6316ac24",
    "metadata": {},
    "outputs": [],
@@ -395,7 +385,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 46,
+   "execution_count": 61,
    "id": "8408d463-557b-481b-afc1-5fbbbd67306d",
    "metadata": {},
    "outputs": [],
@@ -412,7 +402,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 47,
+   "execution_count": 62,
    "id": "445a722a",
    "metadata": {},
    "outputs": [],
@@ -429,7 +419,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 48,
+   "execution_count": 63,
    "id": "05512922-26e5-425f-90a4-0df7c2279ccf",
    "metadata": {},
    "outputs": [
@@ -448,7 +438,7 @@
      "traceback": [
       "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
       "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
-      "Cell \u001b[0;32mIn [48], line 2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mM2E_visualize_interareal_connectivity\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m visualize_interareal_connectivity\n\u001b[0;32m----> 2\u001b[0m visualize_interareal_connectivity(M)\n",
+      "Cell \u001b[0;32mIn [63], line 2\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[38;5;28;01mfrom\u001b[39;00m \u001b[38;5;21;01mM2E_visualize_interareal_connectivity\u001b[39;00m \u001b[38;5;28;01mimport\u001b[39;00m visualize_interareal_connectivity\n\u001b[0;32m----> 2\u001b[0m visualize_interareal_connectivity(M)\n",
       "File \u001b[0;32m~/MAM2EBRAINS/./figures/MAM2EBRAINS/M2E_visualize_interareal_connectivity.py:236\u001b[0m, in \u001b[0;36mvisualize_interareal_connectivity\u001b[0;34m(M)\u001b[0m\n\u001b[1;32m    233\u001b[0m X, Y \u001b[38;5;241m=\u001b[39m np\u001b[38;5;241m.\u001b[39mmeshgrid(x, y)\n\u001b[1;32m    235\u001b[0m ax\u001b[38;5;241m.\u001b[39mset_xticks([i \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m0.5\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m np\u001b[38;5;241m.\u001b[39marange(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;28mlen\u001b[39m(area_list) \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m1\u001b[39m)])\n\u001b[0;32m--> 236\u001b[0m \u001b[43max\u001b[49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mset_xticklabels\u001b[49m\u001b[43m(\u001b[49m\u001b[43marea_list\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrotation\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m90\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43msize\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m6.\u001b[39;49m\u001b[43m)\u001b[49m\n\u001b[1;32m    238\u001b[0m ax\u001b[38;5;241m.\u001b[39mset_yticks([i \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m0.5\u001b[39m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m np\u001b[38;5;241m.\u001b[39marange(\u001b[38;5;241m0\u001b[39m, \u001b[38;5;28mlen\u001b[39m(area_list) \u001b[38;5;241m+\u001b[39m \u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m1\u001b[39m)])\n\u001b[1;32m    239\u001b[0m ax\u001b[38;5;241m.\u001b[39mset_yticklabels(area_list[::\u001b[38;5;241m-\u001b[39m\u001b[38;5;241m1\u001b[39m], size\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m6.\u001b[39m)\n",
       "File \u001b[0;32m/srv/main-spack-instance-2305/spack/opt/spack/linux-ubuntu20.04-x86_64/gcc-10.3.0/py-matplotlib-3.6.2-lhkot3cmeebfk5dp74dnubweq56upksc/lib/python3.8/site-packages/matplotlib/axes/_base.py:73\u001b[0m, in \u001b[0;36m_axis_method_wrapper.__set_name__.<locals>.wrapper\u001b[0;34m(self, *args, **kwargs)\u001b[0m\n\u001b[1;32m     72\u001b[0m \u001b[38;5;28;01mdef\u001b[39;00m \u001b[38;5;21mwrapper\u001b[39m(\u001b[38;5;28mself\u001b[39m, \u001b[38;5;241m*\u001b[39margs, \u001b[38;5;241m*\u001b[39m\u001b[38;5;241m*\u001b[39mkwargs):\n\u001b[0;32m---> 73\u001b[0m     \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[43mget_method\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;28;43mself\u001b[39;49m\u001b[43m)\u001b[49m\u001b[43m(\u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43margs\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n",
       "File \u001b[0;32m/srv/main-spack-instance-2305/spack/opt/spack/linux-ubuntu20.04-x86_64/gcc-10.3.0/py-matplotlib-3.6.2-lhkot3cmeebfk5dp74dnubweq56upksc/lib/python3.8/site-packages/matplotlib/axis.py:1968\u001b[0m, in \u001b[0;36mAxis._set_ticklabels\u001b[0;34m(self, labels, fontdict, minor, **kwargs)\u001b[0m\n\u001b[1;32m   1966\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m fontdict \u001b[38;5;129;01mis\u001b[39;00m \u001b[38;5;129;01mnot\u001b[39;00m \u001b[38;5;28;01mNone\u001b[39;00m:\n\u001b[1;32m   1967\u001b[0m     kwargs\u001b[38;5;241m.\u001b[39mupdate(fontdict)\n\u001b[0;32m-> 1968\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m \u001b[38;5;28;43mself\u001b[39;49m\u001b[38;5;241;43m.\u001b[39;49m\u001b[43mset_ticklabels\u001b[49m\u001b[43m(\u001b[49m\u001b[43mlabels\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mminor\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mminor\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[38;5;241;43m*\u001b[39;49m\u001b[43mkwargs\u001b[49m\u001b[43m)\u001b[49m\n",
@@ -459,8 +449,8 @@
    ],
    "source": [
     "%%capture captured\n",
-    "from M2E_visualize_interareal_connectivity import visualize_interareal_connectivity\n",
-    "visualize_interareal_connectivity(M);"
+    "# from M2E_visualize_interareal_connectivity import visualize_interareal_connectivity\n",
+    "# visualize_interareal_connectivity(M)"
    ]
   },
   {
@@ -585,61 +575,6 @@
     "Go back to [Notebook structure](#toc)"
    ]
   },
-  {
-   "cell_type": "markdown",
-   "id": "57ff902c-d6ce-4f96-9e4f-8e3e7166ab66",
-   "metadata": {
-    "tags": []
-   },
-   "source": [
-    "## S3. Data Loading and Processing <a class=\"anchor\" id=\"section_3\"></a>"
-   ]
-  },
-  {
-   "cell_type": "code",
-   "execution_count": 12,
-   "id": "f5b58845-4d1a-430f-83f4-402fdf918aef",
-   "metadata": {
-    "tags": []
-   },
-   "outputs": [
-    {
-     "name": "stdout",
-     "output_type": "stream",
-     "text": [
-      "loading spikes\n",
-      "Loading data from file\n",
-      "Computing population rates done\n",
-      "Loading data from file\n",
-      "Computing synchrony done\n",
-      "Loading data from file\n",
-      "Computing population LvR done\n",
-      "Loading data from file\n",
-      "Loading data from file\n",
-      "Computing rate time series done\n",
-      "pop_LvR\n",
-      "pop_rates\n",
-      "synchrony\n"
-     ]
-    },
-    {
-     "name": "stderr",
-     "output_type": "stream",
-     "text": [
-      "python3: can't open file './../Schmidt2018_dyn/compute_bold_signal.py': [Errno 2] No such file or directory\n"
-     ]
-    }
-   ],
-   "source": []
-  },
-  {
-   "cell_type": "markdown",
-   "id": "2da9728d-4481-4a15-b810-d125e39cbe4e",
-   "metadata": {},
-   "source": [
-    "Go back to [Notebook structure](#toc)"
-   ]
-  },
   {
    "cell_type": "markdown",
    "id": "bb71c922",
@@ -647,7 +582,7 @@
     "tags": []
    },
    "source": [
-    "## S4. Simulation Results Visualziation <a class=\"anchor\" id=\"section_4\"></a>"
+    "## S3. Simulation Results Visualziation <a class=\"anchor\" id=\"section_3\"></a>"
    ]
   },
   {
@@ -657,7 +592,7 @@
     "tags": []
    },
    "source": [
-    "### 4.1. Instantaneous and mean firing rate across all populations <a class=\"anchor\" id=\"section_4_1\"></a>"
+    "### 3.1. Instantaneous and mean firing rate across all populations <a class=\"anchor\" id=\"section_3_1\"></a>"
    ]
   },
   {
@@ -680,7 +615,7 @@
     }
    ],
    "source": [
-    "from M2E_VISUALIZATION import plot_instan_mean_firing_rate\n",
+    "from M2E_visualize_instantaneous_and_mean_firing_rates import plot_instan_mean_firing_rate\n",
     "plot_instan_mean_firing_rate(M)"
    ]
   },
@@ -689,7 +624,7 @@
    "id": "e91c436e-db94-4cd7-a531-29c032efeeae",
    "metadata": {},
    "source": [
-    "### 4.2 Resting state plots <a class=\"anchor\" id=\"section_4_2\"></a>"
+    "### 3.2 Resting state plots <a class=\"anchor\" id=\"section_3_2\"></a>"
    ]
   },
   {
@@ -722,8 +657,8 @@
     }
    ],
    "source": [
-    "from MAM2EBRAINS_VISUALIZATION import plot_resting_state\n",
-    "plot_resting_state(M, A, label_spikes, data_path, sim_params)"
+    "from M2E_visualize_resting_state import plot_resting_state\n",
+    "plot_resting_state(M, A, label_spikes, data_path)"
    ]
   },
   {
@@ -733,7 +668,7 @@
     "tags": []
    },
    "source": [
-    "### 4.3 Time-averaged population rates <a class=\"anchor\" id=\"section_4_3\"></a>\n",
+    "### 3.3 Time-averaged population rates <a class=\"anchor\" id=\"section_4_3\"></a>\n",
     "Plot overview over time-averaged population rates encoded in colors with areas along x-axis and populations along y-axis."
    ]
   },
@@ -796,11 +731,51 @@
     }
    ],
    "source": [
+    "%%capture captured\n",
     "# 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']\n",
     "# output = {'pdf', 'png', 'eps'}, optional\n",
     "A.show_rates()"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "id": "15bcb6a1-9117-4538-8582-7a89e06167da",
+   "metadata": {},
+   "source": [
+    "0 V1\n",
+    "1 V2\n",
+    "2 VP\n",
+    "3 V3\n",
+    "4 PIP\n",
+    "5 V3A\n",
+    "6 MT\n",
+    "7 V4t\n",
+    "8 V4\n",
+    "9 PO\n",
+    "10 VOT\n",
+    "11 DP\n",
+    "12 MIP\n",
+    "13 MDP\n",
+    "14 MSTd\n",
+    "15 VIP\n",
+    "16 LIP\n",
+    "17 PITv\n",
+    "18 PITd\n",
+    "19 AITv\n",
+    "20 MSTl\n",
+    "21 FST\n",
+    "22 CITv\n",
+    "23 CITd\n",
+    "24 7a\n",
+    "25 STPp\n",
+    "26 STPa\n",
+    "27 FEF\n",
+    "28 46\n",
+    "29 TF\n",
+    "30 TH\n",
+    "31 AITd"
+   ]
+  },
   {
    "cell_type": "markdown",
    "id": "ef74ca3e-98dc-49c9-a4a0-2c640e29b1d9",