{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "b1331599",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Down-scaled multi-area model"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "edec8345-aec1-419e-b9e3-7f612aff8262",
   "metadata": {},
   "source": [
    "![Model overview](model_construction.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f4a649cc-3b68-49e4-b2b6-6f29f13a6d9c",
   "metadata": {},
   "source": [
    "The code in this notebook implements the down-scaled version of spiking network model of macaque visual cortex developed at the Institute of Neuroscience and Medicine (INM-6), Research Center Jülich. The full-scale model has been documented in the following publications:\n",
    "\n",
    "1. Schmidt M, Bakker R, Hilgetag CC, Diesmann M & van Albada SJ\n",
    "   Multi-scale account of the network structure of macaque visual cortex\n",
    "   Brain Structure and Function (2018), 223: 1409 [https://doi.org/10.1007/s00429-017-1554-4](https://doi.org/10.1007/s00429-017-1554-4)\n",
    "\n",
    "2. Schuecker J, Schmidt M, van Albada SJ, Diesmann M & Helias M (2017)\n",
    "   Fundamental Activity Constraints Lead to Specific Interpretations of the Connectome.\n",
    "   PLOS Computational Biology, 13(2): e1005179. [https://doi.org/10.1371/journal.pcbi.1005179](https://doi.org/10.1371/journal.pcbi.1005179)\n",
    "\n",
    "3. Schmidt M, Bakker R, Shen K, Bezgin B, Diesmann M & van Albada SJ (2018)\n",
    "   A multi-scale layer-resolved spiking network model of\n",
    "   resting-state dynamics in macaque cortex. PLOS Computational Biology, 14(9): e1006359. [https://doi.org/10.1371/journal.pcbi.1006359](https://doi.org/10.1371/journal.pcbi.1006359)\n",
    "<br>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b952d0ea",
   "metadata": {
    "tags": []
   },
   "source": [
    "#### Notebook structure <a class=\"anchor\" id=\"toc\"></a>\n",
    "* [S0. Configuration](#section_0)\n",
    "* [S1. Parameterization](#section_1)\n",
    "    * [1.1. Parameters to tune](#section_1_1)\n",
    "    * [1.2. Default parameters](#section_1_2)\n",
    "* [S2. Multi-Area Model Instantiation and Simulation](#section_2)\n",
    "    * [2.1. Instantiate a multi-area model](#section_2_1)\n",
    "    * [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. 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)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d782e527",
   "metadata": {
    "tags": []
   },
   "source": [
    "## S0. Configuration <a class=\"anchor\" id=\"section_0\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "9d6cc7d9-3110-4d96-9f9a-9ec7dee6d145",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create config file\n",
    "with open('config.py', 'w') as fp:\n",
    "    fp.write(\n",
    "'''import os\n",
    "base_path = os.path.abspath(\".\")\n",
    "data_path = os.path.abspath(\"simulations\")\n",
    "jobscript_template = \"python {base_path}/run_simulation.py {label}\"\n",
    "submit_cmd = \"bash -c\"\n",
    "''')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "96517739",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "              -- N E S T --\n",
      "  Copyright (C) 2004 The NEST Initiative\n",
      "\n",
      " Version: 3.5\n",
      " Built: Jul 12 2023 06:25:27\n",
      "\n",
      " This program is provided AS IS and comes with\n",
      " NO WARRANTY. See the file LICENSE for details.\n",
      "\n",
      " Problems or suggestions?\n",
      "   Visit https://www.nest-simulator.org\n",
      "\n",
      " Type 'nest.help()' to find out more about NEST.\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%matplotlib inline\n",
    "import numpy as np\n",
    "import os\n",
    "import nest\n",
    "import json\n",
    "import sys\n",
    "from IPython.display import display, HTML\n",
    "\n",
    "from multiarea_model import MultiAreaModel\n",
    "from multiarea_model import Analysis\n",
    "from config import base_path, data_path\n",
    "\n",
    "sys.path.append('./figures/MAM2EBRAINS')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "7e07b0d0",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "%%capture captured\n",
    "!pip install nested_dict dicthash"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "1d440c07-9b69-4e52-8573-26b13493bc5a",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "<style>\n",
       "table {float:left}\n",
       "</style>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Jupyter notebook display format setting\n",
    "style = \"\"\"\n",
    "<style>\n",
    "table {float:left}\n",
    "</style>\n",
    "\"\"\"\n",
    "display(HTML(style))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27160ba8",
   "metadata": {},
   "source": [
    "Go back to [Notebook structure](#toc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "df83f5ea-1c4b-44d3-9926-01786aa46e14",
   "metadata": {
    "tags": []
   },
   "source": [
    "## S1. Parameterization <a class=\"anchor\" id=\"section_1\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30655817",
   "metadata": {},
   "source": [
    "### 1.1. Parameters to tune <a class=\"anchor\" id=\"section_1_1\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f67c1ba",
   "metadata": {},
   "source": [
    "|Parameter|Default value|Value range/options|Value assigned|Description|\n",
    "|:-------:|:-----------:|:-----------------:|:------------:|:---------:|\n",
    "|scale_down_to|1.       |(0, 1.0]           |0.005         |$^1$       |\n",
    "|cc_weights_factor|1.   |[1.0, 2.5]         |1.            |$^2$       |\n",
    "|areas_simulated|complete_area_list|Sublists of complete_area_list|complete_area_list|$^3$|\n",
    "|replace_non_simulated_areas|None|None, 'hom_poisson_stat', 'het_poisson_stat', 'het_current_nonstat'|'het_poisson_stat'|$^4$ |"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a2161477",
   "metadata": {},
   "source": [
    "1. `scale_down_to` <br>\n",
    "scale_down_to is the down-scaling factor which defines the the ratio of the full scale multi-area model being down-scaled to a model with fewer neurons and indegrees so as to be simulated on machines with lower computational ability and the simulation results can be obtained within relative shorter period of time. <br> Its deafualt value if `1.` meaning full scale simulation. <br> In the pre-set downscale version, it's set as `0.005`, where the numer of neurons and indegrees are both scaled down to 0.5% of its full scale amount, where the model can usually be simulated on a local machine. <br> **Warning**: This will not yield reasonable dynamical results from the network and is only meant to demonstrate the simulation workflow <br> \n",
    "2. `cc_weights_factor` <br>\n",
    "This scaling factor controls the cortico-cortical synaptic strength. <br> By default it's set as `1.0`, where the inter-area synaptic strength is the same as the intra-areal. <br> **Important**: This factor changes the network activity from ground state to metastable state. <br>\n",
    "3. `areas_simulated` <br>\n",
    "This parameter specifies the cortical areas included in the simulation process. Its default value is `complete_area_list` meaning all the areas in the complete_area_list will be actually simulated. <br>\n",
    "complete_area_list = `['V1', 'V2', 'VP', 'V3', 'V3A', 'MT', 'V4t', 'V4', 'VOT', 'MSTd', 'PIP', 'PO', 'DP', 'MIP', 'MDP', 'VIP', 'LIP', 'PITv', 'PITd', 'MSTl', 'CITv', 'CITd', 'FEF', 'TF', 'AITv', 'FST', '7a', 'STPp', 'STPa', '46', 'AITd', 'TH']` <br>\n",
    "The value assigned to simulation_areas can be any sublist of the compete_area_list specifying areas a user want to include in his/her simulation. <br>\n",
    "4. `replace_non_simulated_areas` <br>\n",
    "The paramter replace_non_simulated_areas defines how non-simulated areas will be replaced. <br> It's set as `None` by default when the parameter areas_simulated is set as full_area_list where all areas will be simulated so that no areas need to be replaced. <br> Other options are: `'hom_poisson_stat'`, `'het_poisson_stat'`, and `'het_current_nonstat'`. `'hom_poisson_stat'` is a manually set parameter which can be tuned. When it's set as 'het_poisson_stat' or 'het_current_nonstat', the data to replace the cortico-cortical input is loaded from 'replace_cc_input_source' which is the firing rates of our full scale simulation results. The differenc between 'het_poisson_stat' and 'het_current_nonstat' is that 'het_poisson_stat' is the mean of the time-series firing rate so that it's static, yet 'het_current_nonstat' is time-varying specific current, which is varying by time. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "60265d52",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Downscaling factor\n",
    "# Value range/options: (0, 1.]\n",
    "# Value assgined: 0.005\n",
    "scale_down_to = 0.005 # Change it to 1. for running the fullscale network\n",
    "\n",
    "# Scaling factor for cortico-cortical connections (chi) \n",
    "# Value range/options: [1., 2.5]\n",
    "# Value assgined: 1.0\n",
    "cc_weights_factor = 1.0\n",
    "\n",
    "# Cortical areas included in the simulation\n",
    "# Value range/options: any sublist of complete_ares_list\n",
    "# Value assgined: complete_ares_list\n",
    "areas_simulated = ['V1', 'V2', 'VP', 'V3', 'V3A', 'MT', 'V4t', 'V4', 'VOT', 'MSTd', 'PIP', 'PO', 'DP', 'MIP', 'MDP', 'VIP', 'LIP', 'PITv', 'PITd', 'MSTl', 'CITv', 'CITd', 'FEF', 'TF', 'AITv', 'FST', '7a', 'STPp', 'STPa', '46', 'AITd', 'TH']\n",
    "\n",
    "# Firing rates used to replace the non-simulated areas\n",
    "# Value range/options: None, 'hom_poisson_stat', 'het_poisson_stat', 'het_current_nonstat'\n",
    "# Value assgined: 'het_poisson_stat'\n",
    "replace_non_simulated_areas = 'het_poisson_stat'"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de11b07f",
   "metadata": {},
   "source": [
    "### 1.2. Default parameters <a class=\"anchor\" id=\"section_1_2\"></a>\n",
    "We try our best not to confuse users with too many parameters. However, if you want to change more parameters and explore the model, you can do so by passing a dictionary to the `default_params` argument of the `MultiAreaModel` class."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "6e4bed8d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Connection parameters\n",
    "conn_params = {\n",
    "    'replace_non_simulated_areas': 'het_poisson_stat', # Whether to replace non-simulated areas by Poisson sources with the same global rate, by default: None\n",
    "    'g': -11., # It sets the relative inhibitory synaptic strength, by default: -16.\n",
    "    'K_stable': 'K_stable.npy', # Whether to apply the stabilization method of Schuecker, Schmidt et al. (2017), by default: None\n",
    "    'fac_nu_ext_TH': 1.2, # Increase the external input to 2/3E and 5E in area TH\n",
    "    'fac_nu_ext_5E': 1.125, # Increase the external Poisson indegree onto 5E\n",
    "    'fac_nu_ext_6E': 1.41666667, # Increase the external Poisson indegree onto 6E\n",
    "    'av_indegree_V1': 3950. # Adjust the average indegree in V1 based on monkey data\n",
    "}\n",
    "\n",
    "# Input parameters\n",
    "input_params = {\n",
    "    'rate_ext': 10. # Rate of the Poissonian spike generator (in spikes/s)\n",
    "} \n",
    "\n",
    "# Neuron parameters\n",
    "neuron_params = {\n",
    "    'V0_mean': -150., # Mean for the distribution of initial membrane potentials, by default: -100.\n",
    "    'V0_sd': 50.} # Standard deviation for the distribution of initial membrane potentials, by default: 50.\n",
    "\n",
    "# Network parameters\n",
    "network_params = {\n",
    "    'N_scaling': scale_down_to, # Scaling of population sizes, by default: 1.\n",
    "    'K_scaling': scale_down_to, # Scaling of indegrees, by default: 1.\n",
    "    'fullscale_rates': 'tests/fullscale_rates.json', # Absolute path to the file holding full-scale rates for scaling synaptic weights, by default: None\n",
    "    'input_params': input_params, # Input parameters\n",
    "    'connection_params': conn_params, # Connection parameters\n",
    "    'neuron_params': neuron_params # Neuron parameters\n",
    "} \n",
    "\n",
    "# Simulation parameters\n",
    "sim_params = {\n",
    "    'areas_simulated': areas_simulated,\n",
    "    't_sim': 2000., # Simulated time (in ms), by default: 10.0\n",
    "    # 't_sim': 1500., # Simulated time (in ms), by default: 10.0\n",
    "    'num_processes': 1, # The number of MPI processes, by default: 1\n",
    "    'local_num_threads': 1, # The number of threads per MPI process, by default: 1\n",
    "    'recording_dict': {'record_vm': False},\n",
    "    'rng_seed': 1  # global random seed\n",
    "}\n",
    "\n",
    "# Theory paramters (theory_params)\n",
    "theory_params = {\n",
    "    'dt': 0.1 # The time step of the mean-field theory integration, by default: 0.01\n",
    "} "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1472e9c5",
   "metadata": {},
   "source": [
    "Go back to [Notebook structure](#toc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de4a6703",
   "metadata": {
    "tags": []
   },
   "source": [
    "## S2. Multi-Area Model Instantiation and Simulation <a class=\"anchor\" id=\"section_2\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1fd58841",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 2.1. Instantiate a multi-area model <a class=\"anchor\" id=\"section_2_1\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "ab25f9f8",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Error in library(\"aod\") : there is no package called ‘aod’\n",
      "Execution halted\n"
     ]
    }
   ],
   "source": [
    "%%capture captured\n",
    "M = MultiAreaModel(network_params, \n",
    "                   simulation=True,\n",
    "                   sim_spec=sim_params,\n",
    "                   theory=True,\n",
    "                   theory_spec=theory_params)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "91649c30",
   "metadata": {},
   "source": [
    "### 2.2. Predict firing rates from theory <a class=\"anchor\" id=\"section_2_2\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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"
     ]
    }
   ],
   "source": [
    "# p, r = M.theory.integrate_siegert()\n",
    "# print(\"Mean-field theory predicts an average \"\n",
    "#       \"firing rate of {0:.3f} spikes/s across all populations.\".format(np.mean(r[:, -1])))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2062ddf3",
   "metadata": {},
   "source": [
    "### 2.3. Extract and visualize interareal connectivity <a class=\"anchor\" id=\"section_2_3\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8a7c09e0",
   "metadata": {},
   "source": [
    "The connectivity and neuron numbers are stored in the attributes of the model class. Neuron numbers are stored in `M.N` as a dictionary (and in `M.N_vec` as an array), indegrees in `M.K` as a dictionary (and in `M.K_matrix` as an array). Number of synapses can also be access via `M.synapses` (and in `M.syn_matrix` as an array). <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "6316ac24",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Neuron numbers\n",
    "\n",
    "# Dictionary of neuron numbers\n",
    "# M.N\n",
    "\n",
    "# Array of neuron numbers\n",
    "# M.N_vec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "8408d463-557b-481b-afc1-5fbbbd67306d",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Indegrees\n",
    "\n",
    "# Dictionary of nodes indegrees organized as:\n",
    "# {<source_area>: {<source_pop>: {<target_area>: {<target_pop>: indegree_values}}}}\n",
    "# M.K\n",
    "\n",
    "# Array of nodes indegrees\n",
    "# M.K_matrix.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "445a722a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Synapses\n",
    "\n",
    "# Dictionary of synapses that target neurons receive, it is organized as:\n",
    "# {<source_area>: {<source_pop>: {<target_area>: {<target_pop>: number_of_synapses}}}}\n",
    "# M.synapses\n",
    "\n",
    "# Array of \n",
    "# M.syn_matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "05512922-26e5-425f-90a4-0df7c2279ccf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "6.8556 3.026424152338414\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAg4AAADwCAYAAABovOx5AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjYuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8o6BhiAAAACXBIWXMAAAsTAAALEwEAmpwYAABVr0lEQVR4nO2deZhU1bW339XNKN1MMoOAqFFBBgFxiIoaiTiiSDTOirPR6JcbjYn3Rk1iHJKbqFFxxCEqRokKKg5Eg3rjgAyCgKLMNKPMY3fTXev7Y59qqmuf6joFTVFNrZfnPHSv2nuffaq766yz9m+tLaqKYRiGYRhGFAp29wQMwzAMw6g7mONgGIZhGEZkzHEwDMMwDCMy5jgYhmEYhhEZcxwMwzAMw4iMOQ6GYRiGYUQmreMgIl1FRIMjJiLfi8iLIlKUjQkahmEYhpE7ZBJxmApcAHwGnAdcu0tmZBiGYRhGzpKJ4/A98C9gcvB9vdqfjmEYhmEYuYykqxwpIl2B+UnmpcChqrpyF83LMAzDMIwcJJOIw+fAIOAuoANw3S6ZkWEYhmEYOUsmjsMqVf0XznEAOGUXzMcwDMMwjBwmE51CBxH5KXBM8P2C2p+OYRiGYRi5zI5oHDYAnwJXq+rCXTc1wzAMwzByjbSOg2EYdY8kh78MlxX1EfDfqposdjYMw4iMVY40jD2bqcA1wAfA+cAnItJm907JMIy6jDkOhrFns1RVn1HVS4AngHbA1bt5ToZh1GHMcTCM/OHt4P/eu3UWhmHUacxxMIz8QYL/TdhkGMYOY46DYeQPJwX/T9+tszAMo06TVcdBREaKyEoRmZHidRGRB0VkjohMF5G+2ZyfYeyBdBCRS0XkaeBKYDnw+G6ek2EYuwERaSIiz4rIEyJywY6Ok+2IwzPA4BpePxk4IDiuAkZkYU6GsSdzKPAYcCLwInCkqq7YvVMyDKO2SPVALiKDRWR28CB+a2AeCoxW1SuBM3b0nFnd4VJVPwryy1MxBHhOXXGJz0SkuYi0V9Vl2ZmhYewZqOoCtmsaDMPYc3kGeAh4Lm4QkULgYdz+UiXAFyIyFugEfBU0q9zRE+ba1tgdgcUJ35cENs9xEJGrcFEJmjRp0u+ggw7KygQNY09k8uTJq1S19e6eh2HkCyed1EtXrd5YY5spkxfMBEoTTI+rarWlxhQP5AOAOao6D0BEXsI9mJfgnIcv2YkVh1xzHMKekEIV4MGb9zhA//79ddKkSbtyXoaxRyMiVj7eMLLIqtUb+PTz39bYpmG9yxsCnwBvqOobGQwf9hB+OPAg8JCInApkMl41cs1xKAH2Sfi+E7B0N83FMAzDMHYNCqppVws2Ag2BM0Skqaq+EHH00IdwVd0MXJbBLEPJtXTMscDFQXbFEcB60zcYhmEYewpxMeOsWUuojJVWHW+/PYWe3f+b7gf+hnvvfYPKWClAB6At7t6YiZhxlz6EZzXiICKjgOOAViJSAtwO1AdQ1UeBccApwBxgC7XgGRmGYRhGDvEMTsw4WbUCgMrKGL+48RXGvnUNHTs149gf/pVTTj0Y3D16c9AvEzHjF8ABIrIvsAT4KW6vmloh21kV56V5XYGfZWk6hmEYhhGZkwb30tWrNqV8ffLk+RmJGTXmHIcvJi6k235707VrMwDOHtabN8dOB9gAPK2qb6aquxD2QK6qT4nI9cC7QCEwUlVnZn7F4eSaxsEwjF1EaWkpxx57LGVlZVRUVDBs2DDuvPNOXnnlFYAeIhIDBqhqqNJYRBbg1lwrgQpV7Z+1yRtGDrBq1UY+m3hnytcbFF5SmsnfhbrlCJaWrKJjx6Kq7zt02ItJX5QArAPOrknMmOqBXFXH4aL4tY45DoaRJzRs2JAPPviAoqIitm3bxtFHH83JJ5/MIYccAm55cHWEYY5X1VVRzrf33o21c+dm1WwFW9aFti1bUT/KkAA02v9Az7Z6ZlhSSPiWHC2abfBs6zc09dt17xp5TmHn37tHl8j9d5bVMxeEWMPLeOyKeZXP+9qzNeh2cOT+301d5NkaF4bPv1OvfTzbvC8Xe7Zuffx2AOVzv/FsX63bHC0dWZVYrKymFs1E5HGiZkHEtrlhYxWgsarviVUi7ve3GNgGjMswq2KXYo6DYeQJIkJRUREA27ZtY9u2bYgIBx98MECNn4Y7QufOzfjw3xdXs+019dXQtnP/0jHyuAe+8aFne6bHFZ6tUMIdh7NOfdezvfPujzzbsEnPRp5T2PkvnfRk5P47y98PudyzVWr4jXdXzGvRuQM8W+d/TIzc/9Ti6z1b96aFoW3/NOkBz3ZO8597tpcnPRjaf8GwIzzbvv/8PFI6sqLEAl3CzqNoZRBhaNeQkpJ1Vd8vKVlDu7aNwCUIXFVLJ6w1ci2rwjCMXUhlZSV9+vShTZs2DBo0iMMPPzyT7gq8JyKTgwJshpFnxIjFSlMeBDf6yNGBWAXEKujbpzVz565hwfzVlJeW8eqrX3PK4K4QRDBE5PRdd02ZY46DYeQRhYWFfPnll5SUlDBx4kRmzAjdby4VP1TVvrg9ZX4mIsfumlkaRo6iisYqUh5EuNEHYsZPS0srOfiQx3nu79OoVxjjz/cMZOiwVznsyGc4c8j+HHxgc8jUEckStlRhGHlI8+bNOe6443jnnXfiGoe0qOrS4P+VIvIarqztR7twmoaRY2iVgDEFaZcW4mLGvr1b64fvneWMFaWcdHw7TvrkJ9sbVtR4nt2KOQ6GkScsXryYIUOGUFlZSXl5OZs3b+aRRx5hzZo14HakrQ88LCKDVXVtYl8RaQLMwqWHKdANGFbT+Qq2rqfxzLeq2bb0Dq9h07DhJ55t2rz9Q9ueWHyiZ3u+ly9uXL2lKHxeDcs92+iFzT1bqosrLGji2Spjm0NaZo+LZjzl2UZ2vzJr56/YFl3cGsZbGx/ybGc2vSFy/62Vschtu47+zDdK1P3gdLuAMZwMxJHqlipqJjOxZZYwx8Ew8oQ1a9ZQWVmJqlJQUEAsFqNVq1YMHz4coAkuzbI7MBnoJiIdgCdV9RRc9br2uHTMAuBuVX1n91yJYeweVGNVAsYURBczqiKVaTXJOSmONMfBMPKE3r17M23aNAC2bNnC0UcfjYjw9ddfA8xU1d4i0h6YAFVLE6cEX88TkaXAcVHTMQ1jj0PTRgkyijhI+ohDTmLiSMPII8KyKlasWAEuV5xgb5g2KbpbVoVhaEXqIyMxo0KssuYjR7GIg2HkEfGsinXr1nHWWWftSFbFUhFpA4wXkW9U1cSRRt4gGkNqS7SoilSkXaowjYNhGLuPxYsXc/HFF7N8+XIKCgro3Lkz77zzDs2bN2f9+vUHi8hXwDIg1VLEJyISLzndGsuqMPKOWlyqUI0SVTCNg2EYu4/169dzxx13MHDgQFauXEmXLl04/fTTKS0tBdisqj1FZDTwfXLfIKtCgOOBrcB4oMZwxYr5xdx/8XHVbJcMfC+0bcuRfrbFua2eC2277MrGnq1+0VbPdvpDp4X2H7S1oWd77jY/K+GL4+eE9l9ykZ/tIeJnFbzR95zQ/qdPeTnUnsw/eoVvZjhvY7FnW1nqV1ls1ci/ToADis/0bNe3b+/Z2jQKf7L+2xz/XG9e6v3KcEeX8KyIr9b5GQwLKtZ7trPb+tcJ8Jf9r/Ns/Vr6c/rfkHYAFSkqakYivaAxoxu95PByRE2Y42AYeUIsFuOmm26isrKSWCxGly5dOOCAA9i8eTNAIxH5DlgBtARIkVUxAedAvGhZFUb+UXuCRrfs4acG1wVMHGkYeUKvXr2YOnUq06dP580332TLli0cfvjh9OzZE2Clqh4AvAp0BJdVETgNqOo8YAlQjts22H/ENIw9HY1BRVnqI9MS0XVUHGmOg2HkGZs2beLss8/m/vvvp2nTpowcORKgtYhMxu3Gl+oxyEpOG3mPxCpTHmSSVaGKVJTXeJCje1XYUoVh5AmLFy/mwgsvZPLkyRQVFbF4sduKONA4FAbHMGB5iiFMHGnkN8HNvpYGQ2JpK16aONIwjN1HYWEhxcXFXHHFFfz+97+nX79+DBo0iJtuuglgKdAXJ3r09rjeEXFkm5bruPac16rZ5k4O3xej/nl+GeCVqy4OaQmNn/6rZ5s+9jjP9vyZ4ds6z5zYx7MNm+KXjP7XwPAy0s26Lg21JxNVBJmKAQfMDrX3Km3k2YZO8EtvtF3fMrT/nC1vebZ/Lbvas7Vr7JfxBvi81N+W+9OJQ/wxvw8XVw7t4AtJ79jPf09/88mBof1jIYHyH7b2t1B/f3m4CLJCw7dbj0akTIjIQ+XyckRNmONgGHnCvHnzeOutt+jZsycTJkxgxYoVvPbaayxfvhygK/AN8DUwBUwcaRjJSG1HHCpr3PciZzHHwTDyhKOPPhoNnrYWLFjAscceyw033MDQoUPp3r27Ao2B/sBREFpyOi6OVEwcaeQjqrWXQhmtjkNOYuJIw8gzksWRI0aMAFisqvsA/w/wixo4TBxp5DmKVGxLeZCJmFEVKipqPnIUizgYRh6xbds2zj77bC644AKGDh0KwLPPPguwLmjyCuAvYlMVgUBVV4rIa5g40sg30usSrACUYRh7DosWLWLAgAGUl5dTUlJCYWEhN954Y3z54hAR+RJoR8jnQiCOnAVswH18dsNlYKQkVllA+fqiarYmTbaEtl2/wa8S2PDRv4W2Lb32Zs+2z7RXPdvB/wivnHjXPh08W/uYXw3wkcnNQvv/7a3/hNqT2fTrVqH2tQs87Sn7jJrm2Soqwz+em7Vc59k6aXfP1qqBL0IE+LT/YM828mu/8mLbRuEiwqJG+3q2Vxe29mztG4T3X7TFFy0uXOWLOxsXhgfEi+v7/csr/eyEdo38awKoCJtW+K9lCLWpS1BIn1WRk5jjYBh5wpQpU1ixYkW84BO33HILDRs2ZNy4cRxzzDGVONGjAs9CSnHkRtwS590mjjTyjtrUJcSXKuog5jgYRp5w5plnVokjAYYMGcJ+++3H0UcfDS6b4jBgEfAXCBVHLgWOU9VUm2AZxp5Nbd7sLR3TMIy6xIIFC5g6dSqHH354ovkYYIWqfpeimwLviYgCj6nq47t6noaRa9SWLkGiFYDKScxxMIw8IzmrIoHzgFE1dP2hqi4VkTbAeBH5RlVNHGnkEbWoS1CiRC+ib9OdRcxxMIw8IVXJ6YA2wJXAXBFpq6q3hAxhJaeN/KZWdQmRnBArOQ0gIoOBB3B18Z9U1XuSXm8GPA90Dub3Z1V9OtvzNIw9jVQlp1esWAGwN/CJqh4bRBSqsSMlpyvKGrB8YadqtokLuoW2PfXITz3b2pnhbRve6gdFGj9wnmdb3CF8JeWr/2vv2Q5YsZ9n+69j/DLYAL+I9fNs0+bt79k++8DP1ADYuq1BiNXPqti4sSikXXgGyv5FfgZFqxRZET0GfOnZDijp7NlaNwqff49VR/j9i/3w/bpt4VkRVw/wS4FPn++//+d0XR/a/9uQ619X7mdQ7BcyJ4AYIaWo14Q29alNXYICFWnHsoiDiBQCDwODgBLgCxEZq6qzEpr9DJilqqeLSGtgtoi8oKp1c+Nyw8gRUpWcnjZtGkAF8CK4Og1gJacNIxlBkexGHHKSbFeOHADMUdV5gSPwEpC8O4oCxSIiQBHOF6ybOSuGkUPES05Pnz6d119/nWbNmnHDDTfw7bffAmwGLhORD0XkMHBZFYHTgKrOA+Ilp0uxktNGPqK4m32qI6PKkUBMaz4y2aY7i2R7qaIjkLiwWgIcntTmIWAsbre+YuBcVfXcMhG5CrgKoHNnP8xmGEY4yeLICvcEVQgcgUvJfFlEuql62wiaONLIb1SjLC9EHQwqbakiCmH7nCZ/OJ0EfAmcAOyH+4D6WFU3VOvkUsEeB+jfv//O7JNqGHlBKnFkaWkpQEtgatC0EdAKP6pg4kjDqHl5IbqYMR5xqBkTR+IiDPskfN8JF1lI5DLgnuBpZ46IzAcOAnxFjWEYkUkljjzooIOYO3fuFlXtIyI/AN4HqhV52hFx5KqtjXnqy57VbD//YbjgsPWTMz3bpl/7paEBiu5O/siAyucae7bNt/y/0P599nrEs7044U3PtvXO8JLRjW/36191DW0ZHeWFnervyzVr4kHPcvOfove+NKNz+Wy8eR/PdthepZ7tiY+PDu1/+wK/FPlDB17j2a6f/WjkOf1WHorWMH3EIbMIQXrHISfJtsbhC+AAEdlXRBoAP8UtSySyCPgRgIi0BQ4E5mV1loaxBxIXR37wwQccc8wxVeLIQw89FKCeiMzA6Y4uUVUVkQ4iMi7oniiOnAi8ZeJIIy/RWOojA02CKmiF1niQiWYii2TVcVDVCuB64F1ciduXVXWmiFwjInGX8ffAUSLyFe7J51dW4tYwdp5U4sjCwkJw0ccYbplwKpg40jA8VKEilvrIWByZ5jBxpENVxwHjkmyPJny9FPhxtudlGPlCsjjy2muv5c477/wKp1n4PfC/wPCQriaONIyalxcy0yTUza0qsr5UYRjGbmTbtm2cffbZXHDBBQwdOhSAtm3bAhBkLz2BcyA8Aqc+XufhtVTtDGNPJd3yAhlGHDQmNR4ZjZdFrOS0YeQJixYtYsCAAZSXl1NSUkJhYSE33ngjy5YtA0BEfgn8CecUVCMQR84CNuCCrN2AYTWdr12zDdx68r+q2VqesCC07ZIL+nq2sq37hrZ9vacfDGneYKhnm/uHtaH9L5voPxCuH36IZ2szMnyF9MnuV3q28kq/cuF1GYjzdpZNv+rk2TYuCxd37vPCQs+2+IIunq39c1+G9n/0oKs925Gd53u2ud+3C+3/r6X+z+r5DeM927MHrQvtH8bnqxp6tutTtH3gB9dGHtcjvryQmswiDunrP1lWhWEYu48pU6awYsUKevZ0mQ633HILDRs25D//+Q9AD+B2XMbE/0DKypEbcZHKu00caeQltVjGQSvqZtDfHAfDyBPOPPNMEms6DRkyhP32249rrrmG559/fisu1XIMsAKqliaqxJEishQ4zsTKRt4SLC/UQAbpmAI1j5WzmONgGHnIggULmDp1Kocffjhjx44F2Kaq01yl95Qo8J6IKPBYUITNMPKLihr/RjIqAKWVaSMOVjnSMIzdT2JWRb169bjrrrvAL8QWhmVVGPlN+ohDZsTSOg45qXGomwsshmFkzOLFixk4cCDt2rVjyZIlLF68mLlz5zJjxgyA3iJSjtvOfpqIhCnbPgnqq7zH9pLThpFHBMsLqY6MsirSjJXDyxgWcTCMPCFVyelly5bRrFmzKaraX0RWA+NVdXli3x0pOb1ifVPuf+9H1Wy393gytG2b68o8W4Ojv0ox8s5Vn//2oOs8293fnOLZROpHHlN1207NKYxMzh9Gg3rhkv3K2CbPduAr5Z6t1Wunhvafv+k9z1a0sJtn21r+ReTzh/HX2eElpy9qdIBnG97C//k1abRfaP9YbCd+VumXFzKKEGhINk5dwCIOhpEnpCo53bRp08RmBQQbz1nJacMIIVaQ+sgAjVbHISexiINh5AnxktPgxJHHHnssN9xwQ/zljiKyGFdW+hYIzaqIl5xWrOS0kY+oRBE0RiTScoSJIw3D2P0kl5wOWBIsVfwaVzvn9pCuJo40jAwjCynRSEsVJo40DGP3ESaOBLj55psBeojIdGAg8JMUQ5g40shrVAWtKEh5ZDyeLVUYhpHLpBJHdu/eHWBmEHEYT8gDxY6II1sXb+Laoz+pZtNt4U9Y5WPXe7apJ54Q2raoaLNn2/+1//NshxT9NLT/g/f55Se+GHKiZ1u/9v+F9l947lTPdkbTGzzb74+YHdq/93u+uDCMBw+4ItTeoMAXPXZs4gsOTzjVL+MM8NLoSz3bzQuneLZ5G98K7X9HF/9aLz3MF0J2+6dfhhrgw6NO82x/m7mPZxvY1hdsAnw0+zvPNrSZP6fNpXND+4/qdZlnO/+rZ0LbhqFaSwWgVOqsONIcB8PIE+LiyJ49ezJhwoQqceSUKVNge8ShDLcnRaqS0xNwDsSLJo408g4FrajxZp9ZVkUORxVqwhwHw8gTUokjmzZtiojEIw5vAK+CiSMNw0fSRRwio6SNXuQs5jgYRp6RQhyJiNwGVAAvpOhq4kgjv0kfcchgLKm9sXYSETkTOBVoAzysqjWup5k40jDyiG3btnH22WdzwQUXMHRote2N9wZOAy7QxJ2wEggiEKjqStzW2yaONPKKeJQg1ZHxeDWMFXU8ERkpIitFZEaSfbCIzBaROSJya83z0NdV9UrgUuDcdOe0iINh5AmLFi1iwIABlJeXU1JSQmFhITfeeCO33XYbQJfg6A5MSu4biCNnARtwn5/dgGE1na+wXgXFbVZXs63/pkto24JCX/DXfbAveAQovrWZZ/v+F109WyX9Q/u//9AZnu2hb9p4tnmdFoX2/6TkUM/2Rou/erZ//NCfJ+x8RciNt7b2bNc+7ifCfPDkeaH9n13vv6/X7f1Dz/bQgdeE9n9hzXLPtm2i70P+rN3hof3/Z2qpZ3v14nf9dq8NDu1/W+efe7b9i/1n4Fv38dsBxLRpqD0StRkl2ImxkiIEHwEPAc8lvF4IPAwMAkqAL0RkLFAI3J003PDgYQDgv4N+NWKOg2HkCVOmTGHFihX07NkTgFtuuYWGDRvy97//HdwSRTnwgoj8W1WvSSGO3IiLVN5t4kgjH9lZXYKIjAROa1RYParw/qLN/Pr/viemcOHBTbmpb0uAViKS6Mg/rqqPq+rrwOsi0gL4M/DPpNMMAOao6rzgnC8BQ1T1blxkMXlOAtwDvK2qfopNEuY4GEaecOaZZ5K4CjFkyBD2228/Fi1aRJBRsQn4papOglBx5FLgOFVdtRumbxi7n9qJODyDixBMjmdVVMaUWz7+ntGndqRDk/oMenUhg7sUAWwGEkM8Z4jI6xEiBB2BxQnflwDhISDHDcCJuHTS/VX10ZouwBwHw8hDFixYwNSpUzn88Jo+SzwUeE9EFHhMVf2CCIaxBxMhEyI0QlBtDNWPRKQrbK8cOXnFFroWN6BLk8YAnNmtGePmbQFogNvyvlpdiOQIQXy8BMImGapdCub0IPBgTReWiDkOhpFnpMqqiIBlVRj5jQqxmos2rVLVcHFN2HCBE7JscyUdm9Sv+r79XvWZ8v3WmrpWixAAycuGJUBiVa1OOAekVjDHwTDyhMWLF3PhhRcyefJkioqKqkpOr1mzBuAAoD7wsIgMVtW1IUN8IiIbgUq2l5w2x8HIL2ou2pTRplTxpYpYbPtume4FIQgahBaUSo4QhEQcvgAOEJF9cRvX/RQ4P918omKOg2HUQWbMmMGsWbMoLd2uUL/44otr7JOq5PQzzzwDTvQ4B/gKuBX4VWLfHSk5HasoZPPq5tVsG9eFRzgKQsooN2jml1EG+M8xh3i22VP9D/MnD90W2n/OmuaerVuRr8q/pZuv/geo95dnPNv69Td7tmWXfhjav/Te6Z6tbGULz7b8266h/Sd96GdVxEKC0N+XhXbn3cPberZhX/hlnBcseju0/51Fh3m2f69q5dk+2fJkaP/7D7jOs02e0sez7VtUEdr/5rn+kv7lrW/0bE99Hznyzp8kYluFWC3tjqlQFb1o17ghSzatq/p+yaZK2jZqABEcEREZBRyHWyYpAW5X1adE5HrgXVwmxUhVnVkrE8ccB8Ooc9x5551MmDCBWbNmccopp/D2229z9NFHp3UcUpWcfv7558GldbUGegB7Ab+yktOGUR1FUK3Rccis5HSwNNG7ZRPmbyxn4YZy2u1VnzEL1vHw0V24e+rytOOpamjeraqOA8ZFnUsmWAEow6hjjB49mvfff5927drx9NNPM23aNMrKUjxeJhAvOT19+nRef/11mjVrxg033MCWLVsApqtqQ1VtjUvLRFWXBk4DQVpXvOR0KVZy2shT0uxm2UxEHheR01P1DyIEn5ZVKn1Hz+KFb9dSoIX8vl8nznt/PgPHzua0fVpwQPFekcbbHVjEwTDqGI0bN6agoIB69eqxYcMG2rRpw7x58yL3N3GkYewgCrHYzkUc4hGCXns30XEnHRzY4IQOLTihQ4uEdtHG2x1YxMEw6hj9+/dn3bp1XHnllfTr14++ffsyYEC06s9hJafbtm0LThiJiLQHVob1tZLThrHzEYftA9U8Vi7vnJn1iIOIDAYewAk2nlTVe0LaHAfcj/swW6WqA7M4RcPIaR555BEArrnmGgYPHsyGDRvo1atX2n6XXXYZo0aNokmTJrz7rivxO23aNDZu3AhwcLAz5iRgTHJfEVmIKxC1DYjhlit+V9P5Sssa8fV3+1ez9Tt8cmjbii2NPNtNz5wV2vbWo6Z6tvLyBp5ta4gN4ILr/+7Zlv3pSs/WonhjaP+1czv557/s356t1TPHh/a//9A+nq1RSMntG78dEdr/0Eb7ebYtZWM926t9fhrav/8H/rL3gtCW4Xy/6QvPtvCc6PVAbvruEc+25b/bebZvXzkz8pgNM3gEHtTk2uiNk9D06ZgZRQjSRC8gwyyNbJHViENC/eyTcTXxzxOR7kltmgOPAGeoag/AL8JuGHmMqvL888/zu9/9jq5du9K8eXMmTpyYtl+/fv0oKytj8+bN9OnThz59+jBs2DAeffRRcI7AANwmN/cAiEgHEYnfZQpxQnABGgJvmTjSyEfSbEqVQcSh5mhDEHFYr6pX5ZLTANlfqqiqn62q5cBLwJCkNucDr6rqIqgKixqGEXDdddfx6aefMmrUKACKi4v52c9+lrbf9ddfz/z589l///358ssv+fLLL1mxYgWnn346wLdAf2Crqq6B6uJI3F4Wx6lqb1Xtoap37YprM4ycRoVYZUHKgwxu9Ol22tzZPTF2Jdl2HMLqZ3dMavMDoIWITBCRySISmmMmIleJyCQRmfT99ybwNvKHzz//nIcffphGjVx4v0WLFpSXl+/QWIcccghjx1aFuX9C9WpzicTLTU8WkZwTaxlGNqjtbbUrY4U1HlhWBRCtfnY9oB/wI6Ax8KmIfKaq31br5Op/Pw7Qv3//lDW4DWNPo379+lRWVuLK1cP3339PQcGOPQOMHDmSn//85wAHA8UEqZghWEaFYUD8hp6K6JoEjSSAzMmsimw7DlHqZ5fgBJGbgc0i8hHQGxdKNYy85+c//zlnnXUWK1eu5LbbbmP06NH84Q9/SNtv+PDhjBkzhs2bN1fZSktL2bBhAzinfhjVd+JLJONy02WVhczf0Lya7ZANRaFtx0/yy/sf23ZzSEto222xZ/vkP0d4ts3bwsWRHzx9pmdz+3ZVZ/nalqH9v17miyPP+uUrnm30MQeG9r92chvP9mR/v0rmq4eGF/RaeFFDz3bvOL8a46cpArGbDrncsz0xz78ZfrQlfA+zA4uHerbehb448uXw0/Ny70s82+EHdvFszRtEj6JVZPDoOH6zLzoVqXEzyO2kv9lHvtErkcSROUm2Z11VP1tEGuDqZyfLgccAx4hIPRHZC7cV6NdZnqdh5CSxWIx9992X++67j1//+te0b9+e119/nZ/8JL2G+NJLL+XZZ5+tZrvpppu4/fbbAWbh0jC9O0hSuekf4kT4NZabNow9kfjNPtWR8Xh1VOOQ1YiDqlaE1c8WkWuC1x9V1a9F5B1gOi7t60lVtQ8pwwAKCgr4r//6Lz799FMOOuigjPqOGDGC999/n7KyMjp16sSdd97J8uXLGT58OMAhuL/LKeAyKrBy04bhkeaGntFSRazm8tWZjZdFsl7HIax+tqo+mvT9n4A/ZXNehlFX+PGPf8w///lPhg4dWqVziMKoUaNYsGABp512GjNmOF/8qKOO4qSTTgIXfewPHAVVxZ6qyk2LSLzctGLlpo08RVWorKU6Dm6Tq7SOg2kcDMPYef7yl7+wefNm6tWrR6NGjVBVRCSuVciIESNG8Ne//pVhw4ZNB+4DngJODGlq4kjDIG3EYbeNlU3McTCMOkZQ6TFjwsSRjz76KB9++CG4gmz3kTodM2NxpGHscShU1pKgUVVqbaxsY46DYdRB1q5dy3fffUdpaWmV7dhjj62xz6WXXsrQoUMZNmxYlW2//fbjgQce4Pjjj58FzAe8VIQkceRWYDwRxJGtenThiklPRLqeCyO1Ss3JO9n/xzvZH571LOdmcFFXTPyFZ2va0i/NnIqdXde9KIO2sze+ulPnOmea/16F4edZpObxlQ/s2GQyxG2rbREHcxwMo47x5JNP8sADD1BSUkKfPn347LPPOPLII/nggw9q7BcmjnziiSe48cYbwUUcegLngokjDSMVtRkliKV3HEwcaRjGzvPAAw/wxRdfcMQRR/Dvf/+bb775Jp5SWSNh4kiAyZMnIyKLgE2q+jqYONIwUlFbWRWafotuMHGkYRi1QaNGjarKTZeVlXHQQQcxe/bsnR22JVDTuoKJI428RzWt45DBjV6iRBxyEnMcDKOO0alTJ9atW8eZZ57JoEGDaNGiBR06dNjh8SoqKgBaAP9I1SaIQKCqK0XkNUwcaeQltStoNI2DYRhZ4bXXXgPgjjvu4Pjjj2f9+vUMHjw4bb+wrAqAG264AVxBtndF5C1VvSXx9UAcOQvYgFuq6IYrT10jS79axB1dbqhmO6Jt+CpH/YJKz3Zgt/mhbcNy3xsX++WpF83vHNp/5MwfeLaHl/viuiUX9A3t37Boq2crmeefq8/4d0P7L73oUM+25vu9PVvnl38T2r/gruc8W9G9JaFtw5h/9pGebd9/fhq5/3Vtb/Rsj6yILk4c09dXjR494AvPdt7ffxTa/73Nvmj0pg7+nO5fumsEk7V1s1dqVy+RTcxxMIw6zMCBAyO3XbhwISJSTRzZrVs3xowZA7BYVXsESxGpxJEbcYWi7jZxpJGPaC2mY5J+2SNnMcfBMPKE999/3xNHnnPOOTz33HMMGjToe3BLEcH/yeLIpcBxqrpqN03fMHKA2tMlKEJl+pLTOUndnLVhGLXCt99+y8cffwxwkIh8KCKHpWiqwHsiMllEck7lbRjZQIGYFqQ8Mh6vjm5yZY6DYdQxfvWrX0WyRaGiooK1a9cCfAPcDLws4Rtg/FBV++JqLf1MRGquNmUYeygxlZRHbY6VyxkXtlRhGHWM8ePHc++991azvf32254tmTBxZGlpKS+++CK4AlCPA42AVvi1GjIuOS1Ag0KtZitAQ9uuL2vs2VJtABT2JDbr6wM925byhqH9+7TwxY0i9T3bt2ccEdo/jL0a+2Ouv2nf0LarV+7v2Zava+HZut79TGj/2G2Xeraw+Yukev9inq1evdaerXmjrqH9O2t3z1a/ni/urKgM3zvl/gOu8GzH71Xm2Vo28K8JoF6h/179soNf+7KwsFlof7Qi3B6B2iwTXZdLTtfNWRtGHjJixAh69uzJ7Nmz6dWrV9Wx77770qtXr7T9L730Up59tnq534MOOohDDz0UXNbEOUAFUE3HkFRy+ofAAiKUnDaMPZHaXFqoq0sVFnEwjDrC+eefz8knn8yvf/1r7rnnnip7cXExLVu2TNs/rOT0oYceyltvvQXQA3gJuERV1UpOG4aPQjpBY/TKkenHymi8bGIRB8OoIzRr1oyuXbsyatQoFi9ezAcffECXLl2IxWLMnx9e8yCRUaNGMXHiRHr06EFJSQmXX345hYWFrFmzBtzn2JfAVHBZFYHTgKrOA+Ilp0uxktNG3lKbmoSax8pljYM5DoZRx7jzzju59957ufvuuwEoLy/nwgt3bH/Ja6+9lrlz54JbqlgG/G+KpiaONPIeFyWQlAdByemo0YEISxUZjZctzHEwjDrGa6+9xtixY2nSpAkAHTp0YOPGjWn7DR8+nH79+jFnzpwqW9u2bSksLAS3HPEbwC8r6PhERL4C3mO7ONIw8gtNq3FoJiKPi8jpaYfSmp2QygzHyyamcTCMOkaDBg0QEeJZk8klpFNx6aWXMnToUIYN214tetmyZfG9KpoCa4Gvk/sliSO3AuOJII4sqxTmbCisZjuirV9aGuDYQ/2Sw60e86aSknoh5aHDSkMD9NnqZ1uUXe7bGvzy48jnz4SeEW2puKWTX155/fqbPVvTpn+MPOaKy3xx7eKFfqYFQP8PnvVsnx93mmdrVrQptP8+vV73bIPuP9WzjTxhWmj/l95c69kuaOm/J8/3GBravyIkk+HimSND24aRRpeQ0W6WEZYjcnJ3TIs4GEYd45xzzuHqq69m3bp1PPHEE5x44olceeWVafuNGDGC4cOHV4kjn3rqKW655RYOOeQQgIa4VMz/AVdyWkTGBV0TxZETgbdMHGnkI5pe45BRhECRGo9Mx8sWFnEwjDrGL3/5S8aPH0/Tpk2ZPXs2v/vd7xg0aFDafqNGjfJKTrdu3ZqWLVvy4IMPzsDVb1gBoSWn4+JIxcSRRh5TWUvbaitCRaxuRhzMcTCMOsigQYMiOQs1sWXLFu666y7ee+89HnzwwXTNf6iqS4NNsMaLyDeqattqG3mFKzld480+o/TJIKpQ57ClCsOoYxQXF9O0adNqxz777MNZZ53FvHnzIo8zd+5c5s+fT+/evcEts3cCpohIu+S2QQQivgnWa5g40shH1IkaUx1kkAXhdtqUGg9sqcIwjNrgF7/4BR06dOD8889HVXnppZdYvnw5Bx54IMOHD2fChAmh/ZJLTvfs2ZOrr746vq12JRADTlbV5Yn9AnHkLGAD7qGrGzCMNBQINKqujaRFk3DBXIPiaALPVOx90ELPJg3LQ9tWbmji2QqKtu3U+bPJfSUPRGq3Yd0tofamze/zbG2fnu7bMpjT4RPejNy2/K97eba/HLbCs7VqvzLymAUhD+7nTX86cv+LJZo4MmLRpsjE0kcccnKpwiIOhlHHeOedd7j66qurIg9XXXUV48aN49xzz41vWBXKwoULEZFq4sibb76Z6dOng3MMtuA2ukoljhSgELjbxJFGflJ74kjFSk4bhpElCgoKePnll6vSKkePHl31WvjGlo7333/fE0cm8WegM4SKI5cCx6nqqrCOhpEPxAtA1UAGEQJJNxbkaMlpcxwMo47xwgsvcOONN3LdddchIhxxxBE8//zzbN26lYceeijj8W677TaAXrh0zONTNFPgPRFR4DFVfXxH528YdRl/b9EdI4LQEnJ0qcIcB8OoQ1RWVjJixAjeeCP84ePoo4/OeMy77rqLP/7xj9OBfwLXA7eHNLOsCiPvcYLG2tnkCmrPCck2WXccRGQw8ABurfRJVb0nRbvDgM+Ac1V1dFgbw8g3CgsLmTx58g71TRZHAtx8881xJ6Q7bjvtzoQ7Dp+IyEaciDJecrpGx6H1Xlu57tCZ1WwuYOHT9M+LIl/HkpAqka37lno2Lasf2v+50UM82zXfPBb5/GF8N+QYz9aosT8ngH1e8qtkhhF2nQBlWxt7tm6v/sezhYkgATZs+I3fNoMqk09294uNXTHricj9dWsDz9a+rS+EbNRyQ+Qxm/tD7jLCf4OriF7HIb0TArZUASJSCDwMDAJKgC9EZKyqzgppdy/wbjbnZxh1gUMPPZQzzjiDn/zkJ1X7VQAMHRpeYjdOWMnp7t27c/fdd1O/fv1ZOGfe+yTb0ZLThrGnoUBFLYoW0zghYEsVgHtKmRNs04uIvAQMwSm6E7kBFzY9LLvTM4zcZ82aNey999588MEHVTYRSes4jBgxgvfff78qq+LOO+9k3Lhx/O///i+4iMNMgr9FEemAiwieQvWS0wK8aFkVRn5Se9kOte2EZJNsOw4dgcUJ35cAhyc2EJGOwFnACdTgOIjIVcBVAJ07d671iRpGrvL009Hz0xMJKzl9+eWXAyAis4DlwKtgJacNI4wIWRUZkNsplzWRbcch7F1KjtbcD/xKVStrSi0LVN2PA/Tv3z9CxMcw9gxKS0t56qmnmDlzJqWl29fRR46MvsNfCO1wTv0LKV43caRhALFavNvUVXFktgtAlQD7JHzfCVia1KY/8JKILMBVp3tERM7MyuwMow5w0UUXsXz5ct59910GDhxISUkJxcXFafsNHz6cfv36MWfOnCrbK6+8QseOHcFFA/+sqqk+Fj8Rka+A99gujjSMvELVRRxSHWRYACpCyemcJNsRhy+AA0RkX2AJ8FPg/MQGqrpv/GsReQZ4U1Vfz+IcDSMnqaiooF69esyZM4dXXnmFMWPGcMkll3D++edz0kknpe0fJo5cu3YtjRs3BtgMlIX121Fx5Dcbt3HUh9WfC35SFL4x1y9O/pFn6/H2+6FtO74wJd2pAfhpixtD7SNvGOXZVl8zwbPt/ejsSOcBOGDMx5HbRqXTi1+F2lX98thrrvuBZ9ur/erQ/k3/x7eHlaf+8tTvQvuXVkYrRr31zlah9sa3r/NsXULajT/8nND+Yb9B8zZVerbyPxeF9m/wy/Cy51FJEyXISMwYIXiRk1kVWY04qGoFLk/8XeBr4GVVnSki14jINdmci2HUNQYMcA/59eu7NMPmzZszY8YM1q9fz4IFC9L2HzFiBMOHD69Wcvq+++6jrKwMoDHwgog8CilLTk8AJgJvmTjSyEfiGocaIg7Rx1InjqzpIINNs7JJ1us4qOo4YFyS7dEUbS/NxpwMoy5x1VVXsXbtWv7whz9wxhlnsGnTJn7/+9+n7ZdGHLkFuEBVJ4GJIw0jFSkX83ZorNxdjqgJqxxpGHWElStX8pe//AXYnlnxs5/9DKBaUaddhIkjjbynNlMolborjjTHwTDqCJWVlWzatIkw/WJNGUi1QRCBQFVXishrRKgcaRh7IrUZcdiZ1M5Ae/QRcLuqRt/XvBYwx8Ew6gjt27fnt7/97Q73Dys5vWbNGs4991yAJsDDIjJYVavtzR18QM0CNuAelLrhMp5qpKHsRbeC6mWTV5b6IjaAQf9X4dk2NDogtG0sRBy4TwO/PHOTWBPPBlB4ZjfP1nvgXp6t6fPh4rzFldM829Zyv2Ry00ZdQ/tvKl/m2SpjWz1bmybhZWxaF/n25tLds7Wq3Du0/zd39/Jsh4m/x8noxft6NoAftfm7ZxvbxJfLreLk0P5r/uzbG+KX0b6kdYfQ/n/Y62rPtl78iNsJvz0vtH/st1eE2qMhxEKrCmQwgshI4LR6UlgttfO7rQt5Z83HxFD6FnXnmGb90g31K+DlnZrMDpLtdEzDMHaQ1JmS0Vi4cCEiUk0cOXz4cCZOnBhv0h2YDCnFkYIrS323iSONfMSlY6Y+IvIMMDg+ntuzIsa4NR9yQZvT+Vn785mx+VtWlq8BaCwibyYdbUTkRJwzv6K2rk1EmojIZBE5LV1bizgYRh3h/ffD0xMz6Z8sjrzvvvv45ptv6NChwxTgdFzmRJg4cilwnKqu2qlJGEYdJ00BqFYiMinh+8eTt6BX1Y9EpCu6famipGwlLeo1p2lhcwC67/UDvt4yH2Crqno3chE5Hhcl7A4cGEQxVqrqIQltIm0omUDkCIY5DoZRR2jZsmWtj7lixQrat28PgKouC8SPYSjwnrjtLR9L/jA0jHwgQsnpVaraP+p4cXHk+srNFBcWVX1fVFjE0vLlqeehehuAiFyKiwiOB56Lv55qQ0mCiGHScMOBXrgIRqMo8zbHwTCMKFhWhWGQVhyZUcGm+PJGTLVqGcR9X5WqWWMEQ1WfARCRrklDh24oqap3A+kiGFtFZJyqpkz6MMfBMPKYtm3bsmzZMoA2IvINsJeI3KSq9yc1/YGIfA3Mxz3h/IY0WRVNCupxRFF1gd6hLUpD255Tb3/Ptmhzz9C2y0sLPdvZ+5Z4trELwiM09//EF531FV/u9dApn4f2/8vHgz1bmFgs1adumD60TSP/bpRqzfyMbvM82y0T/WqOTQrCP94f6eZf/9/mlHu2pi0fCe2/2teBckEXv//CWPKmx46Lin0hZrtG/pty8v7fhvb/wxLfvvHFDZ6t3rCNof13BhdxqL2xYkH0oklBMRsqN1V9v6FyE00Km0CGEYwE0m4oWW0u1SMYq2pyGsDEkYaR15xxxhnce++94PafeAF4EDhNRKpSGoKsisbAx8APgXm4zegMI++I1XCQYaXHuDiyXf22rK1Yx9pt66mIVfLNlu/Yr2E3yGDviySibCgZMh99JkpqpzkOhpEnnHfeeRx55JHMnj27Kqvi1ltvjYsu6wPH4dY/PwQuScqq+BtwLFZy2shjImRVpL3Ri8go4NNKKnl0+VNM2zwTpYDjmx7H6NVjGLny7xzQ6ABa1AtPp41IlA0ldxhbqjCMPGHUKH9zJ4CXX36Z7t27VwDnAKW4bIpJqpqYVXEF8E+gAjhaRHqo6szszNwwcoc0WRVpN7lS1fMA2jZoq+e1+qkbE+jSqCuXJNT+SIxg7MA0024ouTNYxMEw8pyDDz4YYDlOmf0OMA3nICQyBeiiqr1x0YfXszhFw8gJ4hqHnYk4VI2VJnqRaQQDl5ZZIiKXp9pQcqffgACLOBiGAe4hoiHQADgGeCzp9b7AGBGZH3zfSkRaWV0HI9/Y2YhDBmNFGi8ewQixextK1hbmOBhGnrNkyRJwOoYDcSLJGfgRh5bAx6p6mogMAEYDq2saN6awMak69DH7fRfattPBczzb8m+7hrZ1pSSq07rrEs/2qy+bhfZ/9fgv/bYvTAltG8b9kVvuGj499gzP9uGWZD8vM340/JD0jQKO28ffHHXCkvaeTe9YE9q/ftE/PNu8Sf7599rLL8MNMOPUjp7tnOH+e6L6YGj/0rtaeLbG/70utK03JmmzKjJKx4xQDDaj8bKFOQ6GYYATR04FyoCZwBwRuQaqtr0fCAwUkWnAVuCnurM1sA2jrqFpb/aRIw5up820zXZU47BLMcfBMPKcjh07glNhd8I5Be+p6ntJzV4DLgy+Xgusz9oEDSNHiHizz2i8NORkxMHEkYaR56xduxagObAv0AFoIiIXJjUzcaRh4G72qQ4yEUfilvNqOsiwLkS2MMfBMPKcf/3rXwBlqvq9qm4DXgWOSmyjqhtUdVPw9Tigvoi0yvpkDWM3UxlLfZDJjT5aVkVOYksVhpHnBFKFvQP9guLq1Vf74BORocDTuJLTe+EiFDWKI9fFNvDG1vHVbC+8FV6GOIymkVuG0+XZn4faY5V+yeq6xJEfjfVsZff4QtDy1eHvYPGfFnu2Fr3nerbV1xwY2n/K5tmebckFff3zjAgV+7NX0z96tjBp5vjDzwntP+hzfwPHrUU/C20bxsrpB4RYv4jUV0ldSjxT4hGHNNhShWEYucc555wDro5DQ9zDRCVwq4hcExdI4sSRBbhStmuAU0wcaeQjlaopDzIsER0h4pCTSxUWcTAMA2CpqvYXkR8Dt6vqd0Bi7uRrwH6q6u2sZxj5gmrt1nGoq663RRwMw0jkp0B4bWo4UkSmicjbItIjm5MyjNxAa/yX2Ug1b5hVW0siuwJzHAzDAEBEGgBnAK+EvGxZFUbeE0/HTHVkSk3LHjuy9JEtbKnCMAyAQuDfuO2zJ4jIcFX9NOH1jEtOtyxsypnFg6rZXux5WWjbH/Wf7NlWLm0b2raoaJNna9WtxLM9OTxcBFm6rijUnsyKy3qF2ts+PT1S/zF9kzNaHUOmPB+p/6QTTgm1bylr6NmO/Y9fVsNv5Si7r9izbV3sV2Pc+1FfBAnh13XskRs9W9MQESTAhnW3eLavTvfPdeQJ0QSLABd12xK5bed/TPSNL4ftQh1ObUl7Iix7QI4WgLKIg2EY4LbgrY/bGKc3bmOcROIlp/sAVwEbSZNVYRh7Gq7kdI0RgoxGS/cvV7GIg2HkORs2bAAoBvYDXlXVcqDcSk4bhk8a7UHk9MnarkKZTcxxMIw8Z968eeA2tXoDt0wxGbgxcBjiWMlpI+9xEYcaXYcMsyrqpudgSxWGkedUVFSAK+o0QlUPBTYDtyY1M3GkYVB7WRDplj0yX/rIHll3HERksIjMFpE5IpL84YSIXCAi04PjExHpne05GkY+0alTp/iXT4rIl8BJODFkIn2BJcHrfyQQR2ZrjoaRC6gqFRpLeWRKhHRMy6oQkULgYWAQbje+L0RkrKom1qGdDwxU1bUicjLwOHB4NudpGPlEu3btwD0AXaGqn4vIHUCTpGZxceRpIjIAGE0acWQ9gb0bVn9qWrQ5eVhHWKZCeE5FOHs36efZmhW0C207qp+voL+8+HrPNvrXi0L7vzPgp56tXbO1nq13t/DVnO+GHOPZ6tff5tn277kytH9FqZ8vsejcAZ5txcrWof3r1zvKs61c39yz/Ti0N5SFlOwu39LYs808+Ueh/T9d2M2znfOZXxak4L7PQ/t/cOQwz9axuDK07a6gtkSLruR02rFyMqsi2xqHAcAcVZ0HICIvAUOAKsdBVT9JaP8ZbqtfwzB2LRXAYyJSAMwDLjNxpGFUR4GKWizNVJnDmRM1kW3HoSOQuMNKCTVHEy4H3g57QUSuwqWF0blz59qan2HkK/FCdpXAOFVdC5g40jCqUXtpkooSM8chEmFVNkLfORE5Huc4HB32uqo+jlvGoH///nXz3TeM3OEbVe0rIm2A8SLyjap+lPB6XBy5SUROwYkjw7YZNIw9FhWokFpaFtG0GRo5S7bFkSW4QjNxOgFLkxuJSC/gSWCIqlqRGcPY9WwDUNWVuOhCtUVzVd2gqpuCr8cB9U0caeQjsSBSEHZkgqYZK5ejEdmOOHwBHCAi+wJLcBvqnJ/YQEQ6A68CF6nqt1men2HkHZs3bwYoCMTLU3DOffLf5VDgaZx4eS+gOWnEkQ0LY+xfvLWabcmWRjs937ZFvrivD7448MtYSGlhYMHanp7tjY0PebbKfzwX2v+sLz/xbFvLfSHlGU1vCO0/dsPfQu3JpCpZHSZOPGfas37/H1wb2v8nR/rz7zP+3UhzApi0urln+0FJ+8hjbjn+VM/W4G+jPVv5L68M7X9sxwc8W8n7h4a2DWN0n0sit01GUSqpPSFmLjsHNZFVx0FVK0TkeuBdXG38kao6M0mE9Vtgb+AREQGoUNX+2ZynYeQTK1asADgIWITLpliuqu+EiCMLcMuNa4BLTRxp5CMxqXF5IYPKkUplTu+BmZqsV44MwpzjkmyPJnx9BXBFtudlGPlKt27dAL4DVgF3Ab+A6n+XuOWL/VT1tKxP0DByBEWpoKJ2BpO0TkjOYiWnDcMAtzwxHLdnRSqODNIxlwK/VNWZWZmZYeQMLheiBjKqu1BXlyqs5LRh5DlvvvkmuCVBf2/r7VjJaSPvUVyUINVBBpUe43qJmv5lMl42sYiDYeQ5H374IbgS0uUEKdMi8ryqJqrz+gJjRGR+8H0rEWmlqqtSjbuuXHizpLoYsnWjnX9WWbHJF/eF0a3YF+EBNG1QFqn/wn/4IkqAreX/idQ/qggyFUOmPB+57aAmvhBy+H6loW0XLtzHs7WX+p6tR5OzQ/vP2PSSZ3tvwLme7eaON4b2v/bwdZ5tzMtDPNvQo8eE9t9ynj/upjFTPNuyi/uE9h+zeGCoPRpKJX6VzwQyizikX6rIycqRFnEwjDznvvvuA/hSVRsAg4FNQHKaQbzkdB9c4bWNpMmqMIw9DbdQkfofGUYcKtL8y2S8bGIRB8PIc4LspfijTz3cA4VayWnDSKY2Iw5p9RIZjpc9zHEwDAOAYOfL/YGHVfVzIHGXISs5beQ9rmhTbaVjRlqqiDxeNjHHwTAMAFS1j4g0B14TkUNUdUbCy1Zy2jBqOeKQZqwMx8se5jgYRp6zePFigB+IyNe4JYuFOK1DouOQsTjSMPY0FKIsL0QcS9NFL8AiDoZh5CLr168HWKqqB4tIa5zjkPwhFRdHniYiA4DRpBFHFojQqLC6/vqEdhtC257d7OeebXNleGnfQvH3ymvb0C/D3L9wv9D+czeUe7YzQ8pDt2kUXrB2WfH1nq0s5t8AzuwUflN4pcS3hW12NLBVeHnuIfvP8WxH7L2/Zzukvd8O4DefHOjZ/nsfPyuj795rQ/tXvtLcs727xC/j/H14UgevTe/t2U49aJZnG3PjEaH96xd42xtx/Ad9PdvXZ/pjAjTyf1UyQKnUGqMEGd3oIzgOFnEwDCP3iLmbXmcRmY4TRi4E5pg40jCScdUXaiDyjV6jLVVYxMEwjNyjV69eALNUtb+IdAU+Aj5X1fEJzUwcaeQ9ilKptVRymkjLHhZxMAwjdxGRIuCfwE2qmrymYOJIw1BQrTHikOEmV2kjDjmJFYAyDANcxch/Ai+o6qvJL6rqBlXdFHw9DqgvIq2yPEfD2K3EIw6pDoIIQbRlBSWmlTUeWAEowzBykcsuuwzgUGCNqp6U/LqIHAeMBeYFps9xDx01iiPXxdYxdkv1ssEvTV8Q2va8jGYcjaHNfMEjwODuMzzb9bP/bxfMIJzrdsGYYTLODb/sHNr2jXfHe7byPxd5ti2L24T2f/9P3q8IN584wbO1e3ZaaP8wxh9+jmdrXC/8afzHF73m2d478SzPdswHvUL7/9f5fsnyJ99KN8M4aTUOGVFXlyos4mAYeU6/fv3ARRyKROTL4DhFRK6JCySBRUD9oF1vTBxp5CVKTCtSHpmN5DI0ajpyFYs4GEaec/3113PDDTd8BRQEe1FUI4g4LFDV07I8NcPIKZS0GoeMRqutmhDZxiIOhmFE4UgRmSYib4tIj909GcPYPdRixEEjRRxM42AYRp3EMioMA4KsitrZqwLSjgU5qnEwx8EwDIAioKOIzAGeVNV74i+o6gYROU5ExgDzgS4ico+q3lrTgEIB9Qr2qmZbdfXBoW1bPfb1zs7f41/l74efq3OKkoZ7GE3/vChy2wa/3OTbUrStPMyXsrY+ya9SWfr7vUP7N/ofX1Pbq/s3frtmG0P7ly1v4dm+37qXZ2va9I+h/eWpkPtwuw9D2/qkXV7IcHfM2hNaZhNbqjCMPKfSlXbuCCwAugPniUj3+Osi0i748mPgKmA58OvsztIwdj+Korot5ZHpaDUte2S69JFNzHEwjDxn8ODB4DImuuFSLucAdyVkVAwDngaOBR7EMiqMPEY1lvLIaByc0LKmI1cxx8Ew8pyrr74aYLWq1lfVTjgNw5JgjwpU9SHgMmAb0Bj4HxNIGvmJolSkPDIeqwYnJFNHJJuYxsEw8pwUwYNkowkkDQOgtm7oSpTliFCxpYgUAL8HmgKTVPXZ2plUNCziYBh5TqdOnaC6Fq4TUG3vYis5bRhQGxEHERkpIiuVbSixqqMytpXyipWUV6ykIrYxLsJMVcJ6CE6XtA0I2ah91yJ7wlJl//79ddKkSbt7GoZRJ6moqKB+/frlwEHAEuAL4HxVnRlvEwgkV6iqisgAYDQuApHyA0REvsdt0W0YuU4XVW2drpFIgUoNgXpl22RVDasAnjCGHAtsAibXK2zp+qlSGVtPYUExUEBlbAOFBUVUxtbPwmUyJTI8ONaq6mMiMlpVh6Wbe21iSxWGkefUq1cPXEnpd4FCYKSqzoyLIwOtwzDgWhGpALYSQSAZ5YPYMOoW+q6yraZIWyMRSXyKfVxVH682gupHwfb1VToG1QrcAoAAikh9YloGsDWsYquIlADlwbdZV1Ga42AYBriQaLUnpbg4Mvj6IeChrM/KMHIIVR1cq+MFNSFcPQepViMiTVbFq8DfROQY4KPanFMUzHEwDMMwjN1BXByplUAs4fsYvj45oZvqFuDyXT29VGRdHCkig0VktojMERGv8pw4Hgxeny4ifbM9R8MwDMPY1WwXR7rvEsWS2URECkTkLhH5m4hckq59Vh0HESkEHgZOJqRCXcDJuDSvA3BV6kZkc46GYRiGkRW0cnu0gRiqFYHeIXpNiHiWhojMSLLX+JCeREZZGtmOOAwA5qjqPFUtB17CTTiRIcBz6vgMaC4i7bM8T8MwDMOodURkFPApgEvJjGsZCnAOQ1woGZlngGrai1QP6SLSU0TeTDraAAcCn6rqL4Br050w2xqHjsDihO9LgMMjtOkILEtsJCJX4SISAGXJ3laO0gpYtbsnEQGbZ+1RF+YI7oPDMIxdjKqeByAi7wCt4tGG6lTZMsrSSKDqIT0410vAEFW9G9jpLI1sOw4SYktWgERpQ/DmPQ4gIpPS5c7mAjbP2qUuzLMuzBHcPGthjMHAA7iUzidV9R4RGYn7oFoJXA+MwT0MdAZKcfns7YDZuL/9vYHVuIp47YH4tomvAWcGfbrjSl9XAuuADcEY84NzFwDFuAeORcFr64AVuCXQAmAN0DLoOwPoH/RvFcyhQTCf6cH53wYuDsYqCPpUBvPsHJz3MFx9i7ZAI9zn1hhcue4FuB1IC3F3hH2AT4L+xcE56uEemlon9J8O/AZ4BGgTXLcCfwbOCd6jBoFtEdABly7bNLBNAa4DxgbjSmD/dWAvx+1RUgBswYWqVwBdgIZB+/8Fzg2uK37PmIHbu2RhcF3xce/APdC1CNoWBv8vDc5fP2i3NGj3WHBd9YNx/w38IHifCwPb2uA9qgjeF3C/Q8NwvxftEs5/H3B+MGb8Zzg/4bxtgjEOU9VJIrIA2Ij7WVaoav8wG7uA2s7QSCLKQ3oiGWVpZHupogT3BxPHq1AXsY1hGDlEDfqlZ6geRv0YOBEYGNR5uA73AX0+0Ac4SFV7A1fgboDXqGof3Af51wljLAX2UdX2wNXAx0G7KcBfcTf5d4BDVbVx0K4fUAZMxt0wW6pqW9zn4HzgQpyDcUYwrxjwk2DcQuAN4JpgjHGB/XbgJ7ibbpw/4p74ZgHf4m7S64GzgjmsAH4btL0kOF8p7sYP7mZ4PzBLVfsElTobBWMdEYz756DtsGDMWcA/gc3AXOCGwPZb3M20BfBfqloQvPYr3M30Y+CRwH4fzuG4GbgAOBj4v+D66gG34W7GzXCfy38Ofi6Dgv63Av8VtD0D56R9iLvhFwKjgrGbBe/XApxzMgPnEDUDugb9TwrGbIa7+ZcCM4FTA1sBLj24FXBe0PYa3M+9AfBsMN+GwXW2BH6K+x1L3kP7+OB97p/GVpeI9ABe9YLqFlW9XFVvUNWH0w2ebcfhC+AAEdlXRBrgfpBjk9qMBS4OsiuOwOWXL0seyDCMnCJUv6SqH+Ge7qtQ1WWqOiX4disu6tAx0DVtCuz1cJ9PKiKdcDeMJ2uagIg0xT0FP7X9VLouocmPcJUxS4PxG4tIPdxNuRR3o/xMVcfjbvSVuJs9OO3VnbgbZQUu+oGq/h2YmjSVEcD3wdef4W6yilvXXhO/ruD1T3E38hUJthKqOyLgnrh/F7RDVVcm9I+/v+fgHIfEa26Gc7LqA88HtjdwN9dCYCBwd2B/AnfDaaCqr6rq7OA9WBDYn1bVClXdCHwH7Bu8F/GoTD1cdCAeqfkrzpHYFlxzfO5xJ7Aj0AS4S1XLAvsMqkd6NuGcijLce980aLc8eK0+8EHQdgLOQVCcc/ps8Po63M/4u+Ca8oFd+wCuqlk9gFNwnvNc4LbAdg3uyQLcL83DwetfAf0jjHlVtq9jB6/d5pln86wLc6yNeeKefJ9M+P4i4KHg6664G8JxuGWIabjQfw/cw0Nl8Lf+NtAT+BLnUGwN2i7HPfUeh7tRrsbdSNbjnkL/HNhm424Sr+NubKUJ4/YARuKWUlbjHIgK3E13fGD7BnfjPQr35Ko4B+BtYENwLT9lu4ItPm7XoF98uWNKcN4luJv0hYF9GS7SsQ0Xtp8ffM59H7RdgAv9Lwuub1sw5xbB98twDkUZblkk8VzfA5MC26zgPBXBuF2C92IuLtryJi6asiBoMxm3bNA1uOYvE2yfBe//guBck3ERlC24qFDidS3H3bAWBOdaEbStwEVT4te1DHfT75zw/Wbclu7J53osaLMg4brWBWP2SLiu+DJLWTCnisBWCtyL+x2Ljzmb4L6S8B5OJvgbCLPl+hH87GYkfF8veD/3xTmJ04AetXa+3X3BdthhR90/cDf2ZMfhb8HXXXGOQ1OgKLCdAszBPa2fn2D7Lvh6H9xa63XBje47nOPwDm5NvUPQfh7uxnkS7sZdGdykmuJ0Ab+Pj4sTqe6He/r6AOcEfBeMf3lw3vuDm80knFPx16B/ZfD6IcENbG3CuF3Z7jh0CNr1C25kH+EehpLtT+O0BZ/joiRbcTfQHjgtwWGB7Tmc8/AN8CDbl1uWJo1ZgXvg6hC0+2XQ/3bgX8AxwHvBe7UtmG8H3E24TWCfDWwOxozbtuCWFzok/Cy34Rw2CbF/kHBdPYP+axOuq1tgm5B0XV2D931N0phxBy1+XdcE/V9Kuq5pOEelAjg+6bqmsN3xa4NzWq4Mvu+QYJ8W/Cw82+7++0rztzcq+N2Jp1LGf5e9h/TaOmx3TMMwaoO0oVFN2GET95TfEXhVVV8MXq/adVNVFwPvA2fj1qX3Bf6B+2B/VFWXBu0L2B6pKMEJwmK4p6xRQN+gXVPgK1WdCxwJzFfVl3Ch7HG4my+qehPuqfdKnBPyXdA/JiI9cDedAmBlfL64iED8GuPXfFzw/9vqiNtX4yIlg3Fh+n2BF4P5tgmuRXARhI3B9QzAPXG/mtC/Ae5GAe4mKcCa4DyX4ESaG3EOwgBV/Ri33LMUF0VYHbRdgdMBNA7OEc/hXxvYtgD/VtWlIlIfF0HZBDwfv67A/hhO7HlwwnVNDObVNLiu+rhlnI+B/yRc15ig/4fBz25bMObjwc9gUsJ1nRX0nxS/LlX9sTpdzMvB3OI6kkKc7mMDzilB3RLPKpwjU/XzCuyvBWN6NnIYVT1PVduran1V7aSqTwX2car6A1XdT1Xvqs1zmuNgGEZtkFa/JCLtAu2S4JYTtgFPBN8jIj/GfdivDtLLTsQ98Q/F3fDOxd1wrhaR4mCXzkLcE+YMVV2Ou+E1wIVqfwTMCtoVAc8Fu3wuAo4IFOQFuHX+r4M5nBTYOuJueKOC/luB03FLDM2BMYG9IDgnuNS54iC75Nrg+qaKSBMR6RO0aRyMOw9YoKptgEG4KMIq3BO0BO2KcOK/r3E33hMCezFuSaE0GPOE4P9JIlIcvFcDg/71gLkisi9O+zEHl3r7etB2bHDshXOYxgU/j2cD29L4dQVz6IiLsiRe11OBbZ/4deEiQ0/hIkrx6/pr0K446br+EthbJ1zXU8H/pQnXpThHpjjhuvqJSHMRKQp+PzYHP493cUspg3HOzCfBz7cJztGbG8y/OMH+Y+C7EFtdSPXPLrs7zJJhSGYwLpw2B7g15HXBhbPm4MQ1fXN0nhcE85uO+4XunWtzTGh3GM7rH5aL72XQ5jjcuuxM4MNcnCdOpPYGLvQ5E7hsN8xxJC4tckaK13fq74dw/dIXuCfkbbgP9CXBHDQ4Twnupjwv+HpFcO6lQbtpwWv3BT/nGcE5ynE323nB+zoz6F+CuyEtYXuq5fzgHM1wKaEzg9fKg7lOxC0/xOeyKfid16DN98H7sjzBvi3oOzP4XoPXtia0iff/DBcVSLSvxEUT1uOesivZrkkoS2i3FfhD8J6UJ9hXJ/TXYIxlgW1jQruyYO7fJ9g2Bdc7D6ep0ODcm4LzLUloGwteK02ybcWlTiaeqzT42cT7V7Jd/7AsoW8smPf/w/0OxPtvC+aTeP71wVyXJPXfjBOrrgi+18D2V9zfYXy+cX1F/PdgW9B+BW4ZKf73uBGXhdEtsMXttRri31OO3T6BDD6UCoM/1G5sF3t0D/ngiof6jgA+z9F5HgW0CL4+OdvzjDLHhHYf4EK5WXccIr6XzYMPls7B921ydJ6/Ae4Nvm6NW8ttkOV5Hgv0JbXjsNv/fuyww47cP+rSUkVdKVeddp6q+omqxsOb8XStnJpjwA24NcKVIa9lgyjzPB+3Tr4IqqWpZZMo81SgOAgDF+Ech+gF6WsBDUmNTCIX/n4Mw8hx6pLjkKoUdaZtdjWZzuFy3FNeNkk7RxHpiBMiPZrFeSUT5b38AdBCRCaIyGQRuThrs9tOlHk+hFtrXYpTet+oqtndAi89ufD3YxhGjpPtktM7Q62Vq97FRJ6DiByPcxyO3qUzCjl1iC15jvcDv1LVykC7tjuIMs96OBX1j3DCsU9F5DNV/XZXTy6BKPM8CafDOAGXEjheRD5W1Q27eG6ZkAt/P4Zh5Dh1yXGoK+WqI81BRHrhKuGdrKqrszS3OFHm2B94KXAaWgGniEiFqr6elRk6ov7MV6nqZmCziHwE9MaJtLJFlHleBtyjqgrMEZH5wEE4YV6ukAt/P4Zh5Dh1aamirpSrjpKW1hmXk31Rlp+MI89RVfdV1a6q2hWXE35dlp2GSPPE5YAfIyL1RGQv3EYuX5NdosxzES4qgoi0xaXEzcvqLNOTC38/hmHkOHUm4qCqFSJyPS4/txAYqaozReSa4PVHcer/eEW6LbinvFyc529xO/A9EjzR77Id2HZijrudKPNU1a/FbU87HZdm9aSqZjXvOuL7+XvgGRH5Crck8CtVzep22yIyCpfS2ErcNrq3E+xKmCt/P3UNEbkNJ9CtxP3+Xa2qn+/eWRnGrkVc5NQwDMPIBBE5Ele86DhVLRORVrgU251a3hGReqq6SzNusnEOY8+lLi1VGIZh5BLtcfqaMgBVXRV3GkTkRyIyVUS+EpGRItIwsC8IHAxEpL+ITAi+vkNEHheR93AVLtuKyGsiMi04jgraXSgiE0XkSxF5TNx25tUQkd+KyBciMiMYM16Zc4KI/FFEPgRuDKoufhhkI70bT70VkSuD/tNE5J/BEqBhVGGOg2EYxo7xHrCPiHwrIo+IyEAAEWkEPAOcq6o9cUvC10YYrx9uK/LzcdUeP1S3B0NfYKaIHIwrq/xDVe2DWx65IGSch1T1MFU9BJdpdFrCa81VdWAw/t9whd364aqKxvczeDXo3xunF7o84vth5Al1RuNgGIaRS6jqJhHph9uH4XjgHyJyK25/hvkJwudngZ/hUpxrYqyqbg2+PgG3tTeqWgmsF5GLcM7FF0EQoTHhxdmOF5FbcHtNtMSVTn4jeO0fwf8H4nb6HB+MVYgrCw1wiIj8AVeVtQin3TGMKsxxMAzD2EGCm/oEYEIgfL0EV68jFRVsj/Q2Snptc5rTCfCsqv46ZQMX7XgE6K+qi0XkjqTzxM8hwExVPTJkmGeAM1V1mohcyvadPg0DsKUKwzCMHUJEDhSRAxJMfXCbNH0DdBWR/QP7Rbgto8HtHNkv+PrsGoZ/n2B5Q0QKRaRpYBsmIm0Ce0sR6ZLUL+4krBK3Y+SwFOPPBloHAk9EpL64bcPB7T65TNzW1mFLIUaeY46DYRjGjlEEPCsis0RkOtAduENVS3GprK8EUYgY20u33wk8ICIf4zQKqbgRt+TwFTAZ6KGqs4D/Bt4LzjceJ9CsQlXXAU/gypq/jqsx4hHsqTIMuFdEpuGiJEcFL/8P8Hkw/jeR3gkjr7B0TMMwDMMwImMRB8MwDMMwImOOg2EYhmEYkTHHwTAMwzCMyJjjYBiGYRhGZMxxMAzDMAwjMuY4GIZhGIYRGXMcDMMwDMOIzP8HpV+sU/gPNtAAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 493.603x217.903 with 4 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# %%capture captured\n",
    "from M2E_visualize_interareal_connectivity import visualize_interareal_connectivity\n",
    "visualize_interareal_connectivity(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bae85d86-157c-47a2-9826-860b410a440e",
   "metadata": {},
   "source": [
    "Full-scale interareal connectivity is from: <br>\n",
    "Schmidt M, Bakker R, Hilgetag CC, Diesmann M & van Albada SJ <br>\n",
    "Multi-scale account of the network structure of macaque visual cortex <br>\n",
    "Brain Structure and Function (2018), 223: 1409 [https://doi.org/10.1007/s00429-017-1554-4](https://doi.org/10.1007/s00429-017-1554-4) <br> \n",
    "Fig 4. D **Area-level connectivity of the model, based on data in a–c, expressed as relative indegrees for each target area**"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e67f37e9-ec8d-4bb1-bd21-45e966f47ab6",
   "metadata": {},
   "source": [
    "Go back to [Notebook structure](#toc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0c1cad59-81d0-4e24-ac33-13c4ca8c6dec",
   "metadata": {},
   "source": [
    "### 2.4. Run a simulation <a class=\"anchor\" id=\"section_2_4\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "15778e9c",
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture captured\n",
    "# run the simulation, depending on the model parameter and downscale ratio, the running time varies largely.\n",
    "M.simulation.simulate()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd6e3232",
   "metadata": {},
   "source": [
    "Go back to [Notebook structure](#toc)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bb71c922",
   "metadata": {
    "tags": []
   },
   "source": [
    "## S3. Simulation Results Visualziation <a class=\"anchor\" id=\"section_3\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c1d7aa61-e85a-4e6a-9e01-e018413a572b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Instantiate an analysis class and load spike data\n",
    "A = Analysis(network=M, \n",
    "             simulation=M.simulation, \n",
    "             data_list=['spikes'],\n",
    "             load_areas=None)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "38ddd973",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.1. Instantaneous and mean firing rate across all populations <a class=\"anchor\" id=\"section_3_1\"></a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bea30fc8",
   "metadata": {},
   "outputs": [],
   "source": [
    "from M2E_visualize_instantaneous_and_mean_firing_rates import plot_instan_mean_firing_rate\n",
    "plot_instan_mean_firing_rate(M)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e91c436e-db94-4cd7-a531-29c032efeeae",
   "metadata": {},
   "source": [
    "### 3.2 Resting state plots <a class=\"anchor\" id=\"section_3_2\"></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aeae56a4",
   "metadata": {},
   "source": [
    "**Fig 5. Resting state of the model with χ =1.9.** (A-C) Raster plot of spiking activity of 3% of the neurons in area V1 (A), V2 (B), and FEF (C). Blue: excitatory neurons, red: inhibitory neurons. (D-F) Spiking statistics across all 32 areas for the respective populations shown as area-averaged box plots. Crosses: medians, boxes: interquartile range (IQR), whiskers extend to the most extremeobservat ions within 1.5×IQR beyond the IQR. (D) Population-averaged firing rates. (E) Average pairwise correlation coefficients of spiking activity. (F) Irregularity measured by revised local variation LvR averaged across neurons. (G) Area-averaged firing rates, shown as raw binned spike histograms with 1ms bin width (gray) and convolved histograms, with aGaussian kernel (black) of optimal width."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ae19bcc3",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "from M2E_visualize_resting_state import plot_resting_state\n",
    "plot_resting_state(M, A, data_path)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "473d0882-8e45-4330-bfa2-2c7e1af0dac4",
   "metadata": {
    "tags": []
   },
   "source": [
    "### 3.3 Time-averaged population rates <a class=\"anchor\" id=\"section_4_3\"></a>\n",
    "An overview over time-averaged population rates encoded in colors with areas along x-axis and populations along y-axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "721d1f03-df25-468d-8075-a807025a9c58",
   "metadata": {},
   "outputs": [],
   "source": [
    "# %%capture captured\n",
    "A.show_rates()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b03d44e8-2216-44ff-ada4-83e9c3e6d30a",
   "metadata": {},
   "source": [
    "|Area index|0|1|2|3|4|5|6|7|8|9|10|11|12|13|14|15|16|17|18|19|20|21|22|23|24|25|26|27|28|29|30|31|\n",
    "|:--------:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|:-:|\n",
    "|Area      |V1|V2|VP|V3|PIP|V3A|MT|V4t|V4|PO|VOT|DP|MIP|MDP|MSTd|VIP|LIP|PITv|PITd|AITv|MSTl|FST|CITv|CITd|7a|STPp|STPa|FEF|46|TF|TH|AITd| "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ef74ca3e-98dc-49c9-a4a0-2c640e29b1d9",
   "metadata": {},
   "source": [
    "Go back to [Notebook structure](#toc)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EBRAINS-23.06",
   "language": "python",
   "name": "ebrains-23.06"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}