Designing excitable behaviors.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Designing excitable behaviors\n",
    "\n",
    "## Overview\n",
    "\n",
    "This notebook provides a guide to designing spiking and bursting behaviors using the simplified neuronal model described in [Neuromodulation of Neuromorphic Circuits](https://arxiv.org/abs/1805.05696) and [Neuromorphic Control](https://arxiv.org/abs/2011.04441). The guide is intended as a general introduction to understanding how spiking and bursting behaviors are built through the interconnection of mixed feedback elements. The elements of this model have a simple form that enables easy synthesis and analysis procedures, but they emulate the same roles that different ionic currents provide in biophysical neurons. By the end of the guide, readers should have a good understanding of the essential elements required for generating a rich set of neuronal behaviors, as well as gain an intuition of what effect different ionic currents have in more complicated biophysical neuronal models.\n",
    "\n",
    "**Note:**\n",
    "When rerunning cells that use `add_conductance` or `add_current` methods, the elements will be added **in addition** to all of the elements added previously. If you would like to change parameters of an element that you have already added, either use element's `update` methods shown later in the notebook, or rerun all of the cells from the start so you start with an empty neuron. You can rerun the whole notebook up to a cell by clicking on the top:\n",
    "\n",
    "`Cell -> Run All Above`\n",
    "\n",
    "or just run the whole notebook from the start:\n",
    "\n",
    "`Cell -> Run All`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining a model\n",
    "\n",
    "We first need to load the model definitions and the required Python packages:\n",
    "\n",
    "- `matplotlib`: Python library for plotting\n",
    "- `numpy`: Python library for numerical operations\n",
    "- `iv_model`: definitions for constructing and simulating neurons"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from iv_model import Neuron\n",
    "\n",
    "%matplotlib inline\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first start by creating a single neuron object of the Neuron class by executing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "neuron = Neuron(v0 = 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This line generates an *empty* neuron circuit, consisting only of a membrane capacitor with default capacitance $C = 1$ with the initial voltage set to $V_0 = 0$.\n",
    "\n",
    "If we instead wish to set the membrane capacitance to say $C = 1\\mu F$ and the initial voltage to $V = 5mV$, we can pass optional keyword arguments when creating the neuron:\n",
    "\n",
    "`neuron = Neuron(C = 1e-6, v0 = 5e-3)`\n",
    "\n",
    "For the purposes of this notebook, we will stick with the dimensionless variables to simplify the exposure, but the parameters of the model may be changed appropriately to generate spiking behavior in the biophysical range."
   ]
  },
  {
   "attachments": {
    "neuromorphic_circuit.png": {
     "image/png": "iVBORw0KGgoAAAANSUhEUgAABiUAAAIHCAYAAADw0F9UAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAuIwAALiMBeKU/dgAAABl0RVh0U29mdHdhcmUAd3d3Lmlua3NjYXBlLm9yZ5vuPBoAACAASURBVHic7N13uCRFuYDx92yCZdkFJOckSRFRQFgUCZIUEMWrIKCACnpREcUsiuDlKoJZuKaLgKACXlQURUAyEgUBSUsGyXlhWWDZc+4fdY47p6e7Z3qmp3rC+3ueepjp7qn6Zg471dNfVxVIkiRJkiRJkiRJkiRJkiRJkiT1i4lVByBJXW5hYC3glcBU4FlguNKIJEmSJEmSJElSX9kIOA2YA4zUlKeA4wmJCkmSJEmSJEmSpJZNBI4mjIYYySkvAZ+sKEZJkiRJkiRJktTjhoBfkZ+MSJZvVhKpJEmSJEmSJEnqaZ+gJuEwbdq0kUMOOWTkkksuGbnjjjtGzjvvvJEDDjhgZPLkycnExDsrjFmSJEmSJEmSJPWYJYAnGU00rLjiiiOzZs0aSXPllVeOTJ8+vTYpcTswucLYJUmSJEmSJElSDzmAmtEPF1xwQWpCYsyJJ56YHC2xY4WxS5IkSZIkSZKkHnIqowmGmTNn5iYkRkZGRubPnz+y8sor1yYljqk0ekmSJEmSJKkHTKg6AEnqEquMPdhss80aHjxhwgQ23XTT1NdLkiRJkiRJSmdSQpKCif9+MHFi3nELXjD+uEklxyNJkiRJkiT1HZMSkhQ8MPbg73//e1MvSBx3f8nxSJIkSZIkSZKkPvVxRteHmDhx4sj111+fu6bEWWedlVzo+h1VBi9JkiRJkiT1gqGqA5CkLrEccBcwFWC99dbjnHPOYaWVVqo78JZbbmGbbbbh4YcfHtv0MLA68EIL7R4DTGnhdT8HrmvhdZIkSZIkSZIkqQscSc3oh2WXXXbk2GOPHXnggQdGhoeHR+66666Rr3/96yMzZsxIjpL4QBttzk3U1Wx5TxttSpIkSZIkSZKkik0G/kKx5MDP2mzTpIQkSZIkSZIkSQNqGnAqjZMCw8D3gUlttmdSQpIkSZIkSZKkAbcTcBUwn/HJgHnA+cDmJbWzCrBaTXkn9QmIUxLHrAYsUlL7kiRJkiRJkiSpS2zP+ATBOh1ubyb1SYkfdbhNSZIkSZIkKYoJVQcgSV3umcTz5yqJQpIkSZIkSeoDJiUkSZIkSZIkSVIUJiUkSZIkSZIkSVIUJiUkSZIkSZIkSVIUJiUkSZIkSZIkSVIUJiUkSZIkSZIkSVIUJiUkSZIkSZIkSVIUJiUkSZIkSZIkSVIUJiUkSZIkSZIkSVIUJiUkSZIkSZIkSVIUJiUkSZIkSZIkSVIUk6oOQJIqtBqwSYNj1ko83xl4ssFrLgYeaTEmSZIkSZIkSZLUh5YFngdGSiyPAFPbiGlmSp0/aqM+SZIkSZIkqWs4fZOkQfYIcHzJdX4HmFtynZIkSZIkSVJfMCkhadAdDcwrqa5ngONKqkuSJEmSJEnqOyYlJA26e4GTS6rr+8DskuqSJEmSJEmS+o5JCUmCo4DhNuuYA3yvhFgkSZIkSZKkvmVSQpLgNuD3bdZxPPBECbFIkiRJkiRJfcukhCQFX2/jtS8R1qaQJEmSJEmSlGNS1QFIUpe4Gjgf2KaF1/4auL/ccEoxCVgdWA5YEZgxuv1R4AbgrjbqXozwWa1KSHA/AJxNWOy7iFWALYDlCcmdWcBfKW/xcUmSJEmSJEmSutL2wEjBMh9Yr8QYZqa08aOCdewKXAY8n1JXbbkTOBiYnFHPoonjAdYFTgFeSKlvLvAFYGKinp1qjrlwdNu2o4+HU+p5cPQ9SJIkSZIkSZLU166hWFKi3bUokspIShyWUkdeOQtYKKWeZFLis4TRDI3q+xkwVFNPbVLiGuDkJuoYBvYr+L4lSZIkSZLU5Zy+SZLGOwo4rcDxR3YqkBI8AfwJuBm4B3gcmAasBLwL2Hr0uLcBnwcOb1DfUaP/fQ74C/BPwuiIVYH3AouP7v8gcAFhREXSRqMF4KHR+O4hTNe0MfAOQt80REjGXALc0fCdSpIkSZIkSZLUgyYAt9HcCIMLO9B+GSMlZgJbUj+NUlLtiIpHU45PjpR4gZC4WCylrqUYP8rkbsJnCeNHSoxNz7QP6Ynx1xGSKWPHntzgPUiSJEmSJEmS1NP2p7mkxI4daLuMpESzJhFGT4y186rE/mRS4jUN6lubMOJh7PixkRi1SYlZhARGnv1qjp8DTG/8ViRJkiRJktQLJjQ+RJIGzkmEu/nzXAucHSGWTnqZMCpkzLINjr+xwf5ZwHk1z3dIOeZBQiIkz8nAU6OPFwG2aHC8JEmSJEmSeoRrSkhSvReBHwL/nXPMd0pqaxrjF4WemnLMJMKohVovEJIKjQwB6wDrAasDSwPLEEYrLMX40Q8LNxdyrotZMIJkgxbrmAdcmajnT23GJUmSJEmSJElS15oOPEn6tE2zaLxeQ7PmZrTRqLynQb0zgK8AdxWo862JOpLTNzXj3TXHXzu6rXb6pgubrOf7Na/5bpOvkSRJkiRJUpdzpIQkpXsWOBY4NGXft4D5ccMp5M3Ar4HlU/bNBe4nTKH0OLAZYeREWWbXPE5bELtZz5ZUjyRJkiRJkrqISQlJyvYD4NOMn9boEeCESqJpzquBP7JgceingJ8D5xLWj7gXGK45/nzKTUrU9itzu6AeSZIkSZIkdRGTEpKU7VHCWga71Ww7kbDmRFlOA6a08Lr7MrZ/kQUJiZsJi03/q4X6W7VkzeMnu6AeSZIkSZIkdRGTEpKU7xTGJyVOKLn+fUqub6eax58kbkICxi+cfVMb9axf8/jmNuqRJEmSJElSF5lQdQCS1OUeSDyfnXpUd1iM8esv/CNy+0PArjXPL26xnpWAjUuoR5IkSZIkSV3GpIQk9Y/hxPMVGhw/EZhRYvt7AmuNPn4K+EOL9XyZEBvAecQf7SFJkiRJkqQOMSkhSf3jWcJC1mP2zTl2E+BSYKMC9S+Zs28b4Nia598Fnks5blFgoYw6hoBPAfuPPh8BjiwQnyRJkiRJkrqca0pIUn85AThs9PEngFcTFtP+F+Gi/2rAu4CtR58X8SBw9mh5EHgeWBl4+2gZq+8q4OsZdWwEPAycDlxOWEx8GFgb2IuQLBnzfeDCgjFKkiRJkiRJktSzNiXcsT9WVqw2nIYWBa5jfMxZ5Y+JY9+aUlcz9dSWS4ElEvXs1EI9/8OCKZwkSZIkSZLUJ5y+SZL6y3PAtsAZOcfcDewN7EJY+6FZBwO3ZOybDRwBbN+gzhuBowgjJNLcC7wPOBCYXyA2SZIkSZIk9YCiU3dI0qDZFLii5vlKwAMVxVLUxsB2hCmWJhKSERcCVxJGIwBMZ8FUfs8B82pevyhhnYoxY33GRoTPZTowB7gVuAR4MSOOnQijMgAuArYCpgBvAjYEJgOPAzcA19TEJkmSJEmSJEnSQOm16ZvKlJy+qVW10zdd2H5YkiRJkiRJ6lVO3yRJkiRJkiRJkqKY1PgQSZLU45YClgSeAR7BKbIkSSrDIsAKhOkvHyZ7KktJktS8ScCyhNkbHiKsYSlJ0kBx+ianb+pVywPHAPcx/u/4GPATYO3qQpMkqed9HJjPgv71JeAvwDtw7UZJklq1OiEJUfsb9ibgi4R1LSVJGggmJUxK9KJ9CAuXj+SUecCheOFEkqRWbEd2H3sesHR1oUmS1LMmsmB0f7I8DGxdXWiSJMUzyEmJRYBrakqr3lRTx09KiEv5Pkd+MiJZflZNmJIk9aTJwOE0Tv7fBSxXUYySJPWitwG30vjmut2qClCSpFgGOSmh3rMd46eSaLZ8pIpgJUnqMQsBf6X5/vVCHJEoSVIzPg4M01z/+hzwymrClCQpDpMS6hUTgBspnpAYW2diRvyQJUnqKT+geB/7rkoilSSpd2wBvEyx/vX0SiKVJCkSkxLqFVvRWkJirHwoesSSJPWONQhTRhTtX8+rIlhJknrIpRTvX+fjNIk9bVLVAWhg/CewYdVBSC1YJvH8KGBOFYFIDWzS5uu/kFLHRwl3rDRrV8I8oM2aA3yqwPEA7wfeWOD4BwlzfxdxEPDqAsffAny3YBtfBlYqcPwVwM8LtnEMML3A8X8Bzihw/ATgfwpFBKcRpj5p1hLANwq28VOKrYOzGuH//yKOAW4vcPxrgQMLtlHUkcB9BY7fhOLJyM8BTxc4flvg3QXb+Ajhh2az/oMwdV2zZgOfKRQRfBB4Q4Hj7wf+q2AbnwLWKXD8Pwl37RfxVWD5AsdfCvyiYBvfIaxJ1ayzgDMLHD8Z+GGhiOCXwEUFjl+a9L/fhrT2+3lrwvpN82u2nQBcXqCOFYDDCrb7PeDmAsevBxxcsI0jgAcKHD8T2LdgG58Gni1w/I7AOwscP0z4rVrE7sA2BY5/Cvh8wTYOADYqcPzdFO8vP0Ox6U/+QfF+/2vU/57KcxHh32wRPwCmFHxNEecCvyn4mh8XPP43o+00azHgmwXb+F/gqgLHrwJ8qWAb3wZuK3D8a4CPFWzjMMJix816E/C+gm0cDMwtcPxOwNsLHP8i4fdGEXsCWxY4/lHC740iil63uwM4umAbXyCce9dalGK/78ZMIHxf1J6T3wR8v2A9hxH62Wb9DTixYBvfBqYVfE0jrXz+nwdWL3D83ym+LujXgVc0OOZ44MqC9UotO4P27uC1WCwWS/yyEMX8V8H6nyhYP4QfVEXa+GcLbfy5YBtFfkSO+UfBNk5qoY1HCrZxVMH6Jxasf4TiF7xWaqGN3Qu2sUkLbby5YBu7ttBG0VLk4hXAe1too8gPNggX2ou2UXT+/aML1v9QwfoBTi7YxrUttFFknYIRwsX8om4u2MZPW2jjqYJtHFGw/qkF6x+heEJwzRbaKFr2LRjT+i20sUPBNrZroY0NCraxTwttLF2wjc8XrH9ewfoh3IhQpI0iCeMxpxZs44oW2ri4YBu/baGN2wu2cVwLbcwp2EbRUvSiHy20cUjB+pdvoY09C7axUQttbF2wjZ1baGPdgm3s30IbixVs4ysF63+uYP0Q/m0UaaPIzTNjfluwjYtbaOPKgm0ULX9uIaaiUyAf30Ibj5f0/tr9/K8o2MapLbRxbxP17g0hqyR10m7Ak4TORpIkSZIkSZI0wJy+SZ02hTD1AoShsedXGIsk9bO1CcOsW/U44e79WsMF67iTYnNnF5mOYczNBdu4t4U2rqPYOVIrd0VfSVhgvFk3tdDGxcDiBY6/o2D9IxSfK73o3aIvtNDGIwWPn91CG08VPP6xFtooanbB4x+meEwvFjz+3hbaKOr2gm082UIbNxVso+i/JQjfI0W+c5Pf1824nGJT7dzSQhsXEqZhaNZdBesfpvj/U/cXPH5uRhtrUmzKAwij9dKmGXmwYD3PZcSUp+hoxCdaaKPo3b4PttDGSwWPv6dgG/MbH1JnVsE2ivT3Y26k8fQXtW5toY1rKPa9fkMLbfyN8DdpVivv4wKKj+4topU7zov+f170fPXFFtooMuURtHZ+VLSPfbyFNopOpfxAC20Umb4WQl9WpI0XCtYP4d9GkTaK9jMQ/o0X6cNvbKGNq6k/Z10U2KxgPWm/XSH8jivqCor9+ygyNeKYi4AZLbwuT6uff5Hf4K20cSmNp+1rZeSyVNgeLBieU/SHrCSpebvR3vDPz8YPWZKknjGTYv3qSxSfekiSpEEzRFgXp0gf+5FKIpXUU0xKSFIcUwl34rWSkHgJWDV+yJIk9ZSLaL5vbWVufEmSBtGHaL5/fYjw21c9zjUlJEnqD3OBr7b42mNpbZojSZIGxfoUmyZxI4ovDi9J0qCZCmxV4PgZwLadCUVSP3GkhCTFMwT8imKjJP5GZ+fhlSSpH7yTsP5As/3rHOD1lUQqSVLvWJqw9kyR37BfrSBOST3GpIQkxbUwcBLNncydS7FFFCVJGmR7As/TuH99GNiiohglSeo1awL/pHH/Op+QkBiqJEpJPcWkhCRVY0fgWtJP5m4H3o8nc5IkFbUKYVTiXOr719nA94DFK4tOkqTeNAX4HCGxn+xfXwb+Cryusugk9RyTEpJUrdcw/oRu92rDkSSpLywCXM2C/vV0YFKlEUmS1PuGgMMYn/BfotKI1BEudC1JUn97OPH8wUqikCSpvzwPPF3z/FHCnZySJKl1I8ADNc/nAU9VFIs6yKSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqSEJEmSJEmSJEmKwqREsCMwseogJEmSJEmSJEnqZyYlYAj4FrBb1YFIkiRJkiRJktTPTErAlsCrgI9VHYgkSZIkSZIkSf3MpAT85+h/twBeXWUgkiRJkiRJkiT1s0FPSqwIvHP08RDw0QpjkSRJkiRJkiSprw16UuLDwOSa5+8HFqsoFkmSJEmSJEmS+togJyUWZsHUTWOmAR+qIBZJkiRJkiRJkvreICcldgeWStn+cWBi5FgkSZIkSZIkSep7g5yU+EjG9lWB7WIGIkmSJEmSJEnSIBjUpMRmoyXLwbECkSRJkiRJkiRpUAxqUuKjDfZvD6wTIxBJkiRJkiRJkgbFICYllgPe0+CYIeDACLFIkiRJkiRJkjQwBjEpsT8wpYnj9gNmdDgWSZIkSZIkSZIGxqAlJSaTvcB10nRg7w7GIkmSJEmSJEnSQBm0pMQ7gRUKHP8xwlROkiRJkiRJkiSpTYOWlPhYwePXA7bpRCCSJEmSJEmSJA2aQUpKvB7YooXXFU1kSJIkSZIkSZKkFIOUlGh2LYmkXYDVywxEkiRJkiRJkqRBNChJiSVpfdHqicCHS4xFkiRJkiRJkqSBNChJif2AqW28/kNtvl6SJEmSJEmSpIE3CEmJibS/LsSSwB4lxCJJkiRJkiRJ0sAahKTE24BVS6jnoBLqkCRJkiRJkiRpYA1CUqLdURJjNgQ2L6kuSZIkSZIkSZIGTr8nJdYFtiuxvrISHJIkSZIkSZIkDZx+T0ocCAyVWN9/ACuWWJ8kSZIkSZIkSQOjn5MSM4B9S65zMrB/yXVKkiRJkiRJkjQQ+jkp8T5gegfqPQCY0oF6JUmSJEmSJEnqa/2alBgCPtqhupcHdutQ3ZIkSZIkSZIk9a1+TUq8BVivg/W74LUkSZIkSZIkSQX1a1LiIx2u/43ABh1uQ5IkSZIkSZKkvtKPSYmVgV0jtONoCUmSJEmSJEmSCujHpMSBwKQI7ewFLBmhHUmSJEmSJEmS+kKMi/cxTQReD/y9wGs2ACaPPr4PeKzAa98CnFbgeEmSJEmSJEmSBla/JSXmAzsUfM1DwHKjj48Cjis1IkmSJEmSJEmSBPTn9E2SJEmSJEmSJKkLmZSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRmJSQJEmSJEmSJElRTKo6AEmSBtRKwLIR2lki8Xwd4PkI7T4N3BmhHUmSJEmS1ENMSkiSVI2DgM9U0O5PI7VzJrBrpLYkSZIkSVKPcPomSZIkSZIkSZIUhUkJSZIkSZIkSZIUhUkJSZIkSZIkSZIUhUkJSZIkSZIkSZIUhUkJSZIkSZIkSZIUxaSqA5AkSXWOBJ6tOogmHQ4sVHUQkiRJkiSpN5iUkCSp+xwLPFR1EE36IiYlJEmSJElSk5y+SZIkSZIkSZIkRWFSQpIkSZIkSZIkRWFSQpIkSZIkSZIkRWFSQpIkSZIkSZIkRWFSQpIkSZIkSZIkRWFSQpIkSZIkSZIkRWFSQpIkSZIkSZIkRWFSQpIkSZIkSZIkRTGp6gAkSZIkSepiCwNrAesASwKLA9NHn495I3AE8CzwDPAwcCtwNzAvZrCSJPWICcDKwNrAaoS+dVFg85pjFgG+DjzNgj72VuB2YHbEWFUykxKSJEmSJAXLA1sBM4F1CcmIVWg8y8BrR0vSPOAeYBZwC3AJcDHh4ookSYNiKiHZ8Gbg1YRExFqExH+ehYHPZ+x7GLiN0MdeDVwA3FFGsOo8kxKSJEmSpEE1FdgZ2IaQjFi35PonEy66rAXsBHwamA9cD1wInAOcN7pNkqR+sinwVmDr0ccLlVz/cqNlS2D/0W3/IiQnLgR+CzxVcptSaR4CRkbLgRXH0o/2YMHn+2LFsUhSN/kmC74fk2X5CuMq6hnS38PvqwxKkqQcQ4Q7NX9GGLGQ1R/HKg8CxwAbdPJNS5IUwarAlwhTLFXdv74A/AZ4OzClk29aaoVJic4yKSFJ6UxKSJIU19LAYYR1Hqq+UJJVrgcOAqZ16DOQJKlsk4D3EkYnzKf6vjStPAb8gDByUeoKJiU6y6SEJKUzKSFJUhwrA98D5lD9RZEiF0++DCzRgc9DkqQyLAR8GLiT6vvNZsvLwK9wdKK6gEmJzjIpIUnpTEpIktRZawH/S/gdUvVFkFbLbOAoYJmSPxtJklo1DTgEeIDq+8lWyzDwB2Czkj8bNcmFriVJkiRJ/WQqYT7rz9DeHNLzCVM93QzcMloeAJ4Dnh397+zRY5cAFgWmEy7WrAqsB7yasHj2Ki3GMB34LHAA8EXgx4QLKZIkVWEXwjRIq7ZZz2PATYS1J24G7iCMaHyWcOPbbEI/PI0F/etiwJLAqwh97HrA2rTW1w8BOwM7AScRzhkea/ndSC1wpERnOVJCktI5UkKSpPK9n9bv3JwLnAt8DtiIcm/imwZsC3wDuBR4qcUYrwdmlhiXJEnNWAc4n9ZHJVxD6AN3IazxVJaJhBsADgBOA55sMcangU+M1idFYVKis0xKSFI6kxKSJJVnacLFiKIXIV4CzgT+g3AXZsx4PwhcRLhYUyTmecB3CCNCJEnqpCHgk4QRDEX72BuATwOrRYx3IcJNACe2GPOlhBGOUseZlOgskxKSlM6khCRJ5dgC+BfFLjrcAxxEuXdrtmo14KvA4xR7D/8gTFshSVInLAWcRbG+aS5wLLBhBfEmTQP2JozSKPIenh19ndRRJiU6y6SEJKUzKSFJUnsmA9+j2EiDqwhTR0yoIN5GphCmn5pFsYs/n6giWElSX9sBeJTm+6MnCNMfzqgi2Ca8ibCwdZFzht8Q1rGQOsKkRGeZlJCkdCYlJElq3QzgTzR/YeE64F10ZzIiaSHgw4TRHM2+v2NxHmxJUjkOJkwV2Ez/8xTwZWCJSiItbgvgLzTfv15P3OmnNEBMSnSWSQlJSmdSQpKk1qwM3ERzFxMeI4w+GKok0vZMJtx1Opfm3ut5xF0XQ5LUXyYQktzN9DnzCaMVeyUZkbQ1cCvNvdeHgY2qCVP9zKREZ5mUkKR0JiUkSSpuDeAOmruIcDqwQjVhlmp9wsKbzbznS+jdC0SSpOpMJiwO3UxfM4uwmHSvmw58m+ZGhTxFGGUhlcakRGeZlJCkdCYlJEkqZgPG/37LKs8Au1cUY6cMEUZNNHPh5EZ661xCklSthYE/0lxC4sfA1GrC7JiZwL00fu/PA2+vKEb1IZMSnWVSQpLSmZSQJKl5awAP0viCwbXAWhXFGMObgPtp/Dn8A1i8ohglSb1jImFB50b9ymzgvRXFGMOShIWwG30OcwlTP0ltMynRWSYlJCmdSQlJkpqzPM1N2fQzwt2e/W4p4Fwafx4X0393s0qSyjME/JTG/cnNwNoVxRjT2KjE+TQekekaE22aVHUAkiRJkiRlWAI4G1izwXFHAId1PhwAXgvsmLHvl4SRDJ30OLAT8HNgz5zjtgBOA3YjTPskSVKtbwAfanDMFcDOwBOdD4cpwCcz9l0HnNPh9keAowhTOZ04Gk+aGcCfCf3sbR2OSX3MkRKd5UgJSUrnSAlJkvJNofECz/OJ/ztu/5x4towYxwTgOzmxjJWfR4xJktQbPknj/uNPwLSIMU3PieW4iHFAWMh7dk48I8A9wLKR4+obE6oOQJIkSZKkFEcDb8zZP58wUiD2hYpuMUy4qPSVBsftS0ikSJIEYX2ibzY45tfArsCczofTlc4DtgGezjlmVeBkvL7eEj80SZIkSVK32Q34eM7+EcKF9lPjhNPVvgYc0+CY7wEbRIhFktTdlgJ+Rf6U/mcB78ep/64BdgGezzlmW+DQOOH0F5MSkiRJkqRusi5hLuehnGM+itMS1foMIfGQZSpwJvCKOOFIkrrQROB0YKWcY/6CaxHVupSQmMibkv4wsteaUgaTEpIkSZKkbjEJOAlYNOeYbwL/EyecnvJp4A85+1cFvh8pFklS9zkE2CpnS0IhpwAAIABJREFU/42EtWFfihJN7zif/GkQJxBulFgmTjj9waSEJEmSJKlbHAJskrP/bOALkWLpNS8T1ti4JeeYvQh3fEqSBss6hOn+sjxJ6B/y1lAYZL8AvpWzfzngO5Fi6QsmJSRJkiRJ3WAF4Es5+x8G9iMs8Kx0zwHvI3+aiW8BC8UJR5LUJY4BpuTs/xhwb6RYetWXCOtMZHkv8MZIsfS8vEVNJEmSVMxkYFfq50G/F7gqfjgt2QhYI7FthDAlSN5FLklq1zHA9Ix9w8DehMSE8v2dMOLkhxn71yJM9XRktIgkSVXaBdg5Z/9PCItfK9+LwO6EfnbxlP1DwLHAxoTRi1Kuhwg/tEeAAyuOpR/twYLP1wsZrVsWeB2wA/BuwqJDWwGvJSxaJ6n3fJMF34/JsnyFcRX1DOnv4fdVBlWhE6j/LO4nDJfuFasBd1L/Pk6uMCZJ/W9rsvvFEbrrAvr+ZMe5ZYVxJf2O7DjnENaYkCT1t6nAXWT3BzcBC1cW3XjTyY7zuArjStqd/HOWj1cXWu9wpITUnZYB3gHsBMwEls45dj5wO3Ad8MfRMrvTAUpSiXYEDmjhdY+3+LpO2Q/YJ7HtQWA74DbgcGDzxP6fAad2IJbpwG8YP1Xnv0ZjbOQeQhL8fGDlmu17AZcCPyonREn6tyHg2zn77yB/Hmyl+zjwFtIXDV8E+C/CVE9lOYowCqOoU+lMXyhJgoOA1TP2jQD/CbwQL5y+cCqwL+F3bJrDgRPx2pwacKREZzlSophXEYbMvUR+1jWvvACcSUhouG6M1L0cKbHAgRl1NCr3l/IOyrEBMJfx8T1HGOU2Zm/q38MVHYrnAyltFb3LeAPg2UQdLwKvLy9MSQLClBJ53/dZP/qr0isjJQA+Q3as84BXltjW1Tlt5ZXDS4xBkrTAIsAjZH//nlBZZOl6ZaQEhP4z+fuvtnyhutDUK0xKdJZJieYsBHyD9pIRaWUW4Q6pvMWMJFXDpMQCvZ6UmArczPjYXqb+Ito06i/yjwBrdyCmv6a0s24L9WxPuGhVW88dpN91K0mtGAKuJfu7/pTqQsvUS0mJieR/vj8vsS2TEpLUXT5F9nfvo8ArqgstVS8lJSAsfJ0V7xNkr5MlvIta6gYrAVcCnyMskJpmGHgMuJWwoM6twJNN1L0W8H3CdFCSpM74CrBeYtt3gbMT2+YAv015/XtLjmclwrpDta4h9B1FnUP9lCprAke0UJckpdma8aPKas0j/OBX6+aTf7fmnoyfqk+S1B+mAAfn7D+K5q4rKdu3CNP1pnkF8MGIsagHOVKisxwpkW8V4G6yM6vnAx8Glst4/STCmhNfIf8OqJU69g4ktcqREgsMEb7Pasv/pdS5buKYiaW8g/a8hjAqojbO68keobYd9e/rdsJnUJZDUto4qI36JlF/9+t8YOP2wpQkICQ/Y9zFX6ZeGikx5iqyY/5uSW1MZHw/vXhKW5dQ3+d7s6QklW8/sr/3H6c7Rz732kgJyB+N8i+ybz4eeC503f3WJswL3atq7xydSJiiqFc9R1gMriyLA38GVkvZdzlh5MQlDep4efTYywl3rb4B+CKwa2lRSlLnjRC+z5LbkuanHFe1rzM+OTJMWHz7pYzj/0o4Oa1NFr8S2IRwwagMeyWezyOsV9SqlwkJ8qtY8F4nEN77dm3UK0lrAttm7HsZF7cu0+HAHzP27UsYTTG3zTbmN3g+ptv6cknqRwfk7PsW4RqX2ncc8GnSbyxckbDe6++iRjRgvkcY8lNbvtqhtpYgzHtW29bVbdTX7SMldqS1eTkt5ZdHGvytijoxo52TaT+T+kbCnbdjdTpSQuo+jpTI95uUOtdqs86yvYX6GE9s4nVHp7zueyXFtH5K3VkXoYr6SUrdbyupbkmD6VCy+8JTK4yrkV4cKTEE3EB23Lt1oM20O14b3XQlSWrfKwk3S6V93z9DuLbajXpxpASEm4Oz4j6twri6WlnDJM8g/A9dWz5EZ6ZV2B1YOtHWmR1oR+qkbYH3p2w/Hngf4a7WdlwGbER7d8ZKkvJ9OfF8Dvnzdo/5Rcq2PShnBGva+hRp7bXiS4QfMbW+WlLdkgbPEOnnw2OS69moPSPAd3L2vy9WIJKkjnsf2dPD/hx4KmIsg+B/gBcy9u1C9yaB+sIQcBf12aCtOtDWJYk2hoFV26jPkRKWZktZIyWGgOtS6r8RmFpSG7X2AhbrQL2S2uNIiXzdPlLi9dTHd3SB11+f8vod2oxpCLgnUedTwMJt1lvrCOrj3rzE+iUNjjeR3Q/eUGFczejFkRIA04BnSY97HrBsye05UkKS4psA3Et2P7V+daE11KsjJQB+SXbs3Xi9uXJljZQYIf2O7OScxu1anTAtTa3LCP/YpF6xHbBhYtswsCftz+Oa5hTq72yVJLXnE4nn84EfFnh92uiFds+bNqf+Ro0zyL5rpxXHUb9exsEl1i9pcOyRs+/X0aIYLHPInmVgEvDOiLFIkjpjJrBKxr4bgX9GjGWQ/DJn3+7RoughZSUlIP3H9X8AC5XYxnupH35U1pQEUiz7p2z7A6FzkCR1v2mEc5xav6XYTRK/pH4B0N2ARduIKy2pUfZ50sPUXyx8O47Ik1Tc9hnbhwlrrKkzTsrZl/U3kST1jrzv8rw+QO35C2EN5DQzCaNAVKPMpMStwDWJbYtT7gKIeyeevwicXmL9UqctQvq/iSJTfkiSqrUj4fu81k8K1vEgcH5i2zTCnKOtmEL9HTj3ABe1WF+e5HtdCNipA+1I6l8rkz0l38XAfRFjGTTnkT0t7VZ0Zl1ISVI822RsHyH/bn61Zx5wasa+ycAWEWPpCWUsqFjrF8DGiW17Ee4ebNeGwHqJbX+g/xdn+SvwiqqDEBC+wNu1KfUXsh4kTEMmSeoNyVESTwEXtFDPLwhT+tXai/QpMRvZgfrzhV9RTt+V9DfCmly1a5+8C3/kSGretjn7zo4WxWCaD5xD+sLWSwAbAVdFjUiSVJYZwGYZ+24gXH9S55wNfDxj37bAnyLG0vXKTkr8CjiGkAEasxNhxMTTbdYdY0qCbjSP/k+8DJJNU7a52Jsk9ZatE8//ALzcQj1nENZoqJ2yaQdgGbKH/maJeZ40QpiT/MM127YmjMAd7lCbkvpL8nu0VitJXhVzAelJCQh32JqUkKTe9Cayr/Xav3bexYTruJNT9r0lcixdr8zpmwAeI8yhVWthwhzJ7ZhIWAS41uPAn9usV4rt1SnbLo8ehSSpVa8Elk1s+0OLdc0BfpfYNgl4d8F6FgN2TWz7O3BLi3E1I/mel6B+RKskZXlzxvZngetiBjKg8qb2y/rbSJK6X94UQRdHi2JwPUf2ecz6hN9MGlV2UgLS78pLu3uviDcDKyS2nUrIPkm9ZJWUbbdGj0KS1Kq04dCXtlFfGedN7yDcBNKo3jJdRv3UUFlDxSWp1gxg1Yx9l+JvvBjuAu7N2PeamIFIkkqV9R0+TGfWmlO9rBEpE0i/UXlgdSIpcSb1UzVtBazURp2DOnWT+s/yKdvandpMkhRPcjTAw6OlVecBDyS2zSR7Adg0yfOkebS2LkURTxMW0q7lSAlJzVg7Z9810aJQ1me9IjAtZiCSpNKsk7H9LuDJmIEMsKtz9mX9fQZSJ5ISLwD/l9LO7i3WN5X6aQxuA65ssT6pSslFrsGOQZJ6yRqJ59e3Wd8w8OuU7c2eN61I/fyk51J8TYpW3JB4nvxsJCnNujn7HEEcT9ZnPYQXTSSpFy0ErJ6xz/41nrzPOu8caOB0IikB5U7h9FbCEN9ap7RYl1S1KSnb5kSPQpLUquSJ/j9LqLOd86Z3U38+F2s06U2J51k/giSpVt4F79uiRaG8z9qkhCT1njUJa/KmsX+N5w5gfsY++9canUpKXAzcndj2OlqbOyv5o3wEOKmVoKQu8GLKtmTSTZLUvZKLXD9SQp3XUz/qYF1goyZemzxPeob6xbM7JTlt1XKR2pXU2/Kmp7s9WhTK+6yLTCEYy5qERPxY6dS1DEnqVXnf3bOiRaEXgfsy9nVj/1qZTnXkI6TPZfzegvUsDuyU2HYZ2YtySd1ubsq2V0SPQpLUquQ822WtC3RyyrZGoyVeBWyc2HYGYSrNGJLv3TnIJTUjbY01CEne2TEDGXB5SYmsv1GVdgBOqymTqg1HkrpO3nf3HdGiEGT3sd3Yv1amk3cXpE0dsDdhjspmvYswJ1qjeqVe8VDKNpMSktQ7kmsDPVNSvadQP8x3D7KHYI/tT4p5npR872nrJklS0vSM7U9FjUJPEW4mTONIbknqPVn9K7iWaWxZn/d0il0X72udTErcClyT2LYqsHmBOpJ3CL4InN5OUFLF7k/Z5sKgktQ7piaeP1dSvQ8CFyS2LQ9snXH8ELBnYtv9wEUlxdOMZxPPJ1J/M4kkJWVdNEl+p6izhsle2y7vwlZVZgN31ZSshIokDaq872772LiyPu8JeCPXv3V6HsZ2Fm5cCdgyse2PeAeNeltyUVCAzaJHIUlqVXJtoGSSoh1Fzps2I8yvXesUwkWmWJLvfRh4KWL7knqTSYnukfWZd2NS4mRCvzdW5lUbjiR1HZMS3SPv83Y04qhOJyV+Rf3Jwu7AlCZeuwf18Tl1k3pdcvQQwPY4J6ok9YrnE8/LPKn8P+pHXryL9Ltp0pIVJ5UYSzMWSzyfi3euSmrMpET3yPrMvWAiSb3HpET3yPu8uzHxX4lOXwh9DDiH8YtVv4JwEfaPDV6b/LH9OPCn8kKTKnEZ4aJN7d2lSwLvAH5TSUSS1J8WBlYYffwi8MDo4yFgJvAmwjnJC8CdwJk0tz7EHML39pgyTyrnAL8jrMFVW//OhEU9x0wG3pN47d+BW0qMpRnJpETWNCCSVCvrBjVHWsWXHP03ppmbCGPImub2JeBfMQORpB6Q991tHxtXVv8KTnf7b50eKQGtTeG0PrBhYttpOERTve954NyU7V/ExW4kqUybEJINdwJ/IHzH7gPMIiSIjwI+BxxGGGHwAPBpGn8XJxctW7y8kIHmzpu2A5Zu4nWdlnzvT1QQg6Tek5XAnBY1CgEsmrG9rPWS2nUHC/ry2vLnKoOSpC6Vd4NQ1ve9OiPv83bUyqgYSYnfA08ntr2d/DsL35uyzamb1C9+mrLtdYQpOiRJ5VuUcD5yAvDKjGOmAUcDhzeo667E8/XaiqzeeSwY1THmrYwfnZFMUswjTJkZ26sTz5OfjSSl6aV1DPqdU2lJUv/I++42KRGXU2k1IUZS4gXgjMS2RQjT1aQZov7H9m3AFSXHJVXlT8CtKdt/BKxcclsTgUOBpUquV5J6yVrALqOPZwFHEKZIej/wbcaPfjgU2DynruSF99eWFOOYYeDXiW2TgXePPl6U+nOoc4FHS46jGcn3blJCUjNMSnSPbk9KfAT48Gj5RsWxSFK3cx2D7mFSogkxkhJQbAqnzYFVE9t+WW44UqWGga+kbF8SOJXyFpZbjrCmy9cIc6tL0iC7lbA2w7qEKZtOIZyfHEK4uH7H6HFDwH/n1JO88L4W5X/H5p03vYP6ha+rGE26EOGzrGVSQlIzsn6MexdnXFPInte6Wy6Y/KSm/L7iWCSp2+VNvWdSIq6sz/slXN/j32IlJS4C7kls25Zw0TQpmawYIcz1LPWT35C+2PtM4BJgpTbqHiJMgXYdsE0b9UhSv7gL2Ag4i3BekfQvwl2YY95M/Q0SY65MPJ8EbNxugAnXAzcktr0RWJ3686RnCItjx7YRYQRHreRnI0lpZmdsXxbXWIsp7bf4mKy/kSSpez2Tsy/vO1/ly/q88/5GAydWUmKE+rmOJwLvSWybAuye2PY36hMaUq8bAfYj/a7SDQgJhc9QbMG/hQgXq64ljC6y05Gk4Bng+QbHnA/cPfp4iJCYSPMP6k8m39l6aJlOTjwfAj5BuKmj1hmEqTJj2zXx/HngmgrikNR77s7YPh1YMWYgAy452q2WI98kqfdk9a+Q/52v8mWtO5j3Nxo4sZIS0NwUTtsDr2jidVI/eBx4G/BIyr6lgG8SEnLHE+YSX4Pxw9oXJyQw9iX8O3mEcBFrw04FLEkdsjjhRLm2pA15XSNxzJolx/G3mscbZBwzDFyV2LZzyXFAmM5vOLHtY4SRGbWqWOAaYKfE86sJC25LUiO35exbJ1oUyvus8/5GeVZlfD+9VsoxU6nv813/TpLal/fdvXa0KLQU9de2x7Tav6oE1xDuEK8ttf8wfpXY9wKwRIdjeqimvQM73JaUZjXgRur/bZRVhoEVYr0ZSU37Jtn/bpevMK6iniH9PRSZ+/nAjDoalftz6tyi5rhrm4zjazWvOaFgvK9uso0izk1pp7bcR9wbTMaslRLLwRXEIak37UL299p/VhhXEfuT/R62rDCuIn5A9ntoNel/dU6deeXwJurerOb4G1uMT5L63WzSv2cvrjKoAqaT3VccV2FcRWxO9ns4tMK4uk7sH7JpC1aPTeE0jXCCWuts4KmORiRV7x5gU+BHwPwS632Z8G/udcCDJdYrSf2qdmHPvMXgTqd+VMAB5YeTet5U69fUj6aIYf/E82HgtArikNSbZuXsc3qJeLI+65eAe2MGIkkqzR0Z2x2JGE/euUzeOdDAiZ2UOJn6H/H7EOZJ3o36+fNPiBCT1A2eJ9wZtj5huqZ25gefBXyJBQuiXt92dJI0GKbUPJ6Tc9xjhFEMtfYlP5HRitOA53L2/7zk9poxFfhgYttFmPyW1Ly7yJ7u7U0xAxlgkwk3RaW5g3BzkySp92RND7QMTuEUS965jNM31UjOS9xpjxJ+xL+tZtsrgY2pX1/iceBPkeKSusWthIs9nyGssbIN8HrCPOppU5nNI/ywu4owF/rfWDAVlCT1iuuAY1p43dMlx1E7p/UTDY79IePPZ2YA7weOLTGeOcCnSJ9GYzZwS4ltNWtP6udI/UEFcUjqXfOAG4CNUvZtSPiOeTJqRINnE7IT6X9vo95fABe28LrL2mhTkrTANcAeGfu2wTv1Y9gmY/scqvn9php7UT+n1i8IJ6e1234cKR7XlFCvGAKWJCQo1iAsJBc7sSipPK4p0VmtrClxUc1r9mni+GsZ/54fISQn+tUihHU8at/zPwn9kyQV8XWy+8B3VRhXs3p9TYnDyI7//RXGlcc1JSSpsdeS/f1+eoVxNavX15RYm+z4vfE+oYoLmv9HuItwsZpte6ccd0KUaLrf6ixYd0PVep5q7wYdIdy52+juXUlSccsTFiUbc1ETr/kOcFLN82WATxAWzO5HHwFWSmz7Lo7Ok1TcRcDnM/ZtRfjNqM7JS5xcGCsISVLpbiSMNkyObIbw3T+E5+6dZP9aQBVJiReAM4D9co6ZBVweJ5yutw7wjaqDEBCmH3OKCknqT19jwXnRucA9TbzmFML0ShvWbPss8L/03xoLSwGHJrbdhDeRSGrNpYSR8pNT9u1CSPAOR40omAm8uYnj3pCzb0/CXf2N/I5q5pZenPFJ+Fq3A/dFjEWSVK5hwsXv3VL2LU34/q9iyrzFgQ83cdxCOfteB3yuiTpuI/SxVdgpZ9/50aJQrq3JHs4yQhhOGku3T9+0I/mflSVeeaTB30qSinD6ps6qnb7pLsKUd2kmMX4akfkUW2h1JuHkv/a9n0f/TWl0JvV/460rjUhSr7uU7H5wy4piWouwyHOnf1c8QfaaDp12QE5cP6oopmY4fZMkNeejdOf3/JU5cZVZqpoGckngxYyYngQmVhSXEiYQ7sBI+0MNk76gY6eYlLA0W0xKSCqTSYnOqk1KjJ1fXAT8F3AQoc8/Grg5cdyRLbR1KvXv/0Pthd9V0tYDq+ruI0n943Nk94M/qzCu03LiKqscEe3d1MtLBu1cYVyNmJSQpOasRv1NU2PlSfJHI3TSbhkxlVluJVxzrkJeMujkimJShs2Ad6eUXSLHYVLC0mwxKSGpTCYlOiuZlGimHE1rIxyWJkzxV1vXbGD9tt5Bd1gHeIrx7+0Jeuv/UUndaQWyRyXMBhapKK7XZ8RUVplDmBKvCuvmxPUg3X0Xp0kJSWreX8n+vn93RTFNoP6GsLLLB6K9m3rX5MS1bYVxqYuZlLA0W0xKSCqTSYnOqk1K3ESYruJG0mO9AnhLm+3tnFLvA8DKbdZbpWUJa2sk39euFcYkqb+cS3ZfuE+Fcf0lJ652y/civo+ko3Li+m6FcTXjzSyI1aSEJOX7ANnf92dXGNd+OXG1W+4HpsR7K+NsmBPXA3R30l8VMilhabaYlJBUJpMSnVWblLi2Zvt6hIvq7yYkEtYqsc0fUf85XAnMKLGNWKaRPsXH8VUGJanvvI/svvAmqpuCodEaiK2Wl4BVIr6PWq8gjEDJiu31FcXVrJ1YEOs/Ko5FkrrdDMLIvLTv+2HCotFVmEL2dP7tlk9GfB9Jp+fEdXSFcXW1SVUHoIYuJu4aG2XbmQV3A80jDBnuVfOrDkCS1LZbRksnfAxYm/ELQL8BuBDYjjDtUS+YQbiDamZi+8XAh+OHI6mPnQEcByyasu9VhCTyb6NGFFxASCpvWnK9vyRcjKnCQWQvrn0T4xP43WjpmsdPVxaFJPWG2YQ14PZM2TcEfAF4T9SIgpeAbwPfKbneJ4Cfllxns9YjrJeR5aRYgfQakxLd73ngrqqDaMOjNY9H6O33IklSnpcJJ/4Hp+zbDDgrbjgt24yQgLg4sf2HhBsMJKksc4ATCEndNJ+nmqQEhKmOziixvhHCKMkqTAc+nrP/2FiBtGH7msedurlAkvrJcaQnJSBcRF8bmBUvnH/7KXAosGSJdf4QeK7E+or4AtkjOy/GKQeVo9unb+p1e7Dg832x4lgkqZs4fVNnZU3fJEnqLqsQ7pzM6hN3qSiuIeCGnLiKlt/EDX+cQ3PiegBYuLrQ6uxAmJt77ALPFOCjhOlGxmLeqZrQJKnnXED29/8vK4wrr18qWmYDS8QN/9/WJdy0lRXb9tkvlUxKdJpJCUlKZ1Kis0xKSFLvOJ7sPvFOYGpFce2dE1fRsknk2MesQva84iNUOwd3mrH/F+YQ/vbJ84zLqG6tEUnqNduR/f0/DGxVUVyLEabiK6N/PSZy7LXOyYnrqgrjUo8wKdFZJiUkKZ1Jic4yKSFJvWMN8u80PLKiuCYBd+fE1Ww5N3bgNc7MietBumuUBOQnqK4AlqkuNEnqSZeS/b06C1ioorjyfg83W14AVowd+Kg9G8TmqD41ZFKis0xKSFI6kxKdtR5w2mg5quJYJEmN/YLsfnEusGZFcX0sJ65my7bRow52aBDXZyuKK8/bgVMId5jeSViE+3fAB4DJFcYlSb3qbeT3BYdUFNdyhP69nf61qsWtpwH35cR1LWEaSCmXSYnOMikhSelMSkiStMDyZPcpI8D5wMQK4lqY8b8Zi5ar44cMhKkx7syJ61aquztWkhRX3qi5ZwmLXlfhJzlxNSrzqS7uH+TENQxsXlFc6jEmJTrLpIQkpTMpIUnSeAeRfwHiqxXFdViDuPLKHhXEC3B6g7jeUlFckqT4ViN/faF/UM10fusRkgut9K+/rSBegHc1iOtnFcWlHmRSorNMSkhSOpMSkiSNNwG4nPy7D99eQVytLsh5E9UsyvypBnGdWEFMkqRqHUJ+3/DziuI6rUFcWecDG1YQ63qEkSVZcT0CLFFBXOpRJiU6y6SEJKUzKSFJUr1Nyb9r8iFghQriamVBzg9WEOcbCAt/ZsX0NL11niFJKscUQrI8r9/au4K4XktIMhTpX8+pIM6pwDUN4qqi31cPMynRWSYlJCmdSQlJktL9N/k/+m8AFo8c07IUW5DzXuIvzLwW4S7NvLj2ihyTJKl7bEy4NpfVR7wAbFNBXH/JiSmtxJ6CcBL563KMAH/Axa1VkEmJzjIpIUnpTEpIkpRuEnAZ+T/+LybctRjTjxvEVFsOiRzbcuQvbD0CHB85JklS9zmY/L7iGeD1kWPapkFMteXyyLENEdaJyIvpfmDJyHGpD5iU6CyTEpKUzqSEJEnZVgGeJP8iwO+AiRFjWgN4uUFMI8BjwLSIcc0ArmsQ082RY5IkdachGt/1/zDwyshx/a1BTGNlt8hxHdkgnpeBLSLHpD5hUqKzTEpIUjqTEpIk5XsHjeeZ/iVhnuxYTm0QzwhweMR4lgAubRDP88BrIsYkSepuSwL3kd933E2YFjCWdzaIZ4SwJkbMKZIObSKmQyPGoz5jUqKzTEpIUjqTEpIkNXYIjS8InA0sGimeTRrE8gJh/YkYVgRubBDPPGCXSPFIknrH+jQekfgIsFGkeCYCdzSI54BIsUwAvt8glhHgBFxHQm0wKdFZJiUkKZ1JCUmSmnM0jS8MXAksFSme83LiODZSDOvw/+3dd7gkVZn48e+dPMQBBoackaSIBAMoAiKCBDFjRMGw6mJOqGvYlRWzLLsoJjCCYEAFQQQVkCCggEiOkgYYYIY4w6T7++Pt/t261af7VnWo7nvn+3meemZudfep07HqnPec98AdLeoxTMwyObyi+kiSxp/didl0rc4ljwL7VFSft7eox73A9ArqMA04pUU96tuZwNQK6qMJzKBEbxmUkKQ0gxKSJBUzBPyAsTsIbqWaEZ0vbnL8JcS6E712MPBQkzpkt6MqqIskaXw7iDh/tTqfPAW8j97PCpgO3NOkDh/p8bEBNgIuanL87HYJrtOkLjAo0VsGJSQpzaCEJEnFTQV+x9gdBYuAf6+gPqlOi5N6fMypwFcYe52NYSLthCRJRRxOsXPLL4FZPa7L+xPHnUfvgwAHAA8mjp3frifW5JA6ZlCitwxKSFKaQQlJksqZCZzO2B0Gw8Bp9Dad0yG54y0Dtunh8bYkRmYWee5fxRzXkqRy3gYsZexzzG3Abj2sx0pEECJ7zM/28HgzibZ5kaDMlcC6PayLVjAGJXrLoIQkpRmUkCSpvCGi071I5/zjREfGtB7V45/cM7tzAAAgAElEQVSZY53eg2MArAYcy9ipNeqBkff0qB6SpInvRcQaEkXOsb8l0h31wmczx3kMWLNHx3kzzdNF5bffYcomdZlBid4yKCFJaQYlJElqzxDw3xTrRBgGLgd27UE93pI5xrN7UP5BxIjUIs9xCXBYD+ogSVqx7AHMp9i552HgXcCkLtdhTSIYMQx8o8tlA2xGtFeLXkf8nGoW2dYKxqBEbxmUkKQ0gxKSJHXmcOAJinUoLANOAXbs4vGnAncAf+5imRAjVf9I8c6SubXHSJLUDdsxejbgWNs/gNcBk7tYh68T/YgbdrHMjYHjgCcpfu3webr7vKT/z6BEbxmUkKQ0gxKSJHVuO6IzpGjHyXLgDGD3Lh3//cD+XShnEvAy4K8Ufy7DwO+BOV04viRJWTOBEyh3TroFeDvdmVWwIfDtLpQDsDVwIrCY4s/lXgz4q8cMSvSWQQlJSjMoIUlSd8wEjqdcx8kwcA3wUWCDDo69Ep0tKr0NMQrzjjHqmt8WAx+j+ykzJEnKejWwgHLnqIeI8/LzOjx2J2s4rE4s3n0BxRaxzm6/A9bp4NhSIQYlesughCSlGZSQJKm7Dgb+RfngxFLgHCId1MYV1HM74H2UnxVR366kN2tYSJKUsilwFu2ds24CPk2s7dTrNEizgVcBPwMWtlHXBcCRdDbYQCrMoERvGZSQpDSDEpIkdd/KwBeItkc7nSfDwK3A94A3A88iZmK0a3XgOcRCoCczuv1ZdptPdJaY21qS1A8vp73gf317BPgt8CFgT2DtDuoyFdiWmMlxLJHKseyMiPq2HPghsG4H9VFJU/pdAUmSJEmSuuQJ4CjgJCJ1xN5tlLF5bTu89vdyIr3S9cA9wOPAY7V/H63dZw1gFWBVIjCyCdFZ0q0FOn8MfBi4v0vlSZJU1q+ItYw+RZyTppZ8/GrAgbWt7iHgOmI9iieI8+sjxPl1GXFOrZ9fVwfWImYcbtnG8VOuIwYOXNCFsqRSnCnRW86UkKQ0Z0pIktR7+xEdDe2O6uzntgz4JZHyQpKkQbI58E1gEf0/X7azXQccRncCG2qDi2JJkiRJkiaqs4E9gBcwkg970C0FfgTsALwCuLy/1ZEkqcFtxAyDzYCvErMHx4O/E+tOPB34AbCkv9XRisyZEr3lTAlJSnOmhCRJ1dueOAffQ/9Haea3G4BPAhv17NlLktQbaxD9qpfQ//NpflsAfIcYoCANDIMSvWVQQpLSDEpIktQ/k4F9iRkJ8+lfR8m9wHHEYtiSJE0ETwP+C7iJ/p1fnyAW1X4tMKO3T1ftcKFrSZIkSdKKZhlwTm2DyI29T23bl1hMsxfuBc4AzgX+QgySkyRpIrkJ+I/atjawJyPn1017dMyngAuJ8+u5wFXEuV4DyqCEJEmSJGlFdxvw7dq2ErAzsA2wFTHic2sicDGtYHlPAjcTHTP1f68nclkv7WbFJUkaYPOA02rbEJFGcXtGzq318+ysguUtA/7FyLn1xtq/VxAzHzVOGJSQJEmSJGnEk8Roywtz+6cCqwArE8GJnzKSdulU4ChgIbAIeBRHaEqSlDUM/LO25a3ByHn2TcBna/sfJQYKLKn9fxFxrtU4Z1BCkiRJkqSxLSFGYdZHYj6Wue1BYraFJEkqr35ufQC4J7N/KXBL9dVRr03qdwUkSZIkSZIkSdKKwaCEJEmSJEmSJEmqhEEJSZIkSZIkSZJUCYMSkiRJkiRJkiSpEgYlJEmSJEmSJElSJQxKSJIkSZIkSZKkShiUkCRJkiRJkiRJlTAoIUmSJEmSJEmSKjGl3xWQJEkNLgGW9bsSBa3S7wpIkiRJkqTxw6CEJEmDZ5N+V0CSJEmSJKkXTN8kSZIkSZIkSZIqYVBCkiRJkiRJkiRVwqCEJEmSJEmSJEmqhEEJSZIkSZIkSZJUCRe6liSpP+4B/lbBcaYAz8z8fSPweAXHvbWCY0iSJEmSpHHGoIQkSf1xbG3rtbWBBzJ/vx24sILjSpIkSZIkNTB9kyRJkiRJkiRJqoRBCUmSJEmSJEmSVAmDEpIkSZIkSZIkqRIGJSRJkiRJkiRJUiUMSkiSJEmSJEmSpEoYlJAkSZIkSZIkSZUwKCFJkiRJkiRJkiphUEKSJEmSJEmSJFXCoIQkSZIkSZIkSaqEQQlJkiRJkiRJklQJgxKSJEmSJEmSJKkSBiUkSZIkSZIkSVIlDEpIkiRJkiRJkqRKGJSQJEmSJEmSJEmVMCghSZIkSZIkSZIqYVBCkiRJkiRJkiRVwqCEJEmSJEmSJEmqhEEJSZIkSZIkSZJUCYMSkiRJkiRJkiSpEgYlJEmSJEmSJElSJQxKSJIkSZIkSZKkShiUkCRJkiRJkiRJlTAoIUmSJEmSJEmSKmFQQpIkSZIkSZIkVcKghCRJkiRJkiRJqoRBCUmSJEmSJEmSVAmDEpIkSZIkSZIkqRIGJSRJkiRJkiRJUiUMSkiSJEmSJEmSpEoYlJAkSZIkSZIkSZUwKCFJkiRJkiRJkiphUEKSJEmSJEmSJFXCoIQkSZIkSZIkSaqEQQlJkiRJkiRJklQJgxKSJEmSJEmSJKkSBiUkSZIkSZIkSVIlDEpIkiRJkiRJkqRKGJSQJEmSJEmSJEmVMCghSZIkSZIkSZIqYVBCkiRJkiRJkiRVwqCEJEmSJEmSJEmqhEEJSZIkSZIkSZJUCYMSkiRJkiRJkiSpEgYlJEmSJEmSJElSJQxKSJIkSZIkSZKkShiUkCRJkiRJkiRJlTAoIUmSJEmSJEmSKjGl3xUYAEcDK9f+f0k/KyJJkiRJkiRJ0kRmUAL+t98VkCRJkiRJkiRpRWD6JkmSJEmSJEmSVAmDEpIkSZIkSZIkqRIGJSRJkiRJkiRJUiUMSkiSJEmSJEmSpEoYlJAkSZIkSZIkSZUwKCFJkiRJkiRJkiphUEKSJEmSJEmSJFXCoIQkSZIkSZIkSaqEQQlJkiRJkiRJklQJgxKSJEmSJEmSJKkSBiUkSZIkSZIkSVIlDEpIkiRJkiRJkqRKGJSQJEmSJEmSJEmVMCghSZIkSZIkSZIqYVBCkiRJkiRJkiRVwqCEJEmSJEmSJEmqhEEJSZIkSZIkSZJUCYMSkiRJkiRJkiSpEgYlJEmSJEmSJElSJQxKSJIkSZIkSZKkShiUkCRJkiRJkiRJlTAoIUmSJEmSJEmSKmFQQpIkSZIkSZIkVcKghCRJkiRJkiRJqoRBCUmSJEmSJEmSVAmDEpIkSZIkSZIkqRIGJSRJkiRJkiRJUiUMSkiSJEmSJEmSpEoYlJAkSZIkSZIkSZUwKCFJkiRJkiRJkiphUEKSJEmSJEmSJFXCoIQkSZIkSZIkSaqEQQlJkiRJkiRJklQJgxKSJEmSJEmSJKkSBiUkSZIkSZIkSVIlDEpIkiRJkiRJkqRKGJSQJEmSJEmSJEmVMCghSZIkSZIkSZIqYVBCkiRJkiRJkiRVwqCEJEmSJEmSJEmqhEEJSZIkSZIkSZJUCYMSkiRJkiRJkiSpEgYlJEmSJEmSJElSJQxKSJIkSZIkSZKkShiUkCRJkiRJkiRJlTAoIUmSJEmSJEmSKmFQQpIkSZIkSZIkVcKghCRJkiRJkiRJqoRBCUmSJEmSJEmSVAmDEpIkTVxbAx/M7XsHsEsf6iJJ0kSyB7B95u89gVcC0/tSG0mSJoa1gUMzf68CvA/YoD/VkTReHQoM17an+lwXSVpRzAFOBpYz8huc384BntavCkqSNM69gfT59U7gtX2slyRJ49kqwGIaz69PAV8HZvavapLGE4MSklStZxAdIs2CEdltPjGyU5IkFbcFcBKtz7HH9KtykiSNUzOBDwMLaX5+vQKY3a8KSho/DEpIUnXWoXhAIhuY2LYflZUkaRw6DFhEsXNsPoWiJElK2wq4gWLn1wuBaf2ppqTxwqCEJFXnRMoFJLIXdZIkqbUDgWUUP78+BWzel5pKkjR+rAncRrk27If7UlNJ44ZBCUmqxibAUtoLSgxjGidJklqZCdxF+fPrCf2orCRJ48hxlD+/PoizJca1Kf2ugFYoQzhSSJJ65c3A5A4e/xYi9VPW7cQFX1FrArNK3H85cEeJ+wOsDaxa4v6LgbtLHmNdYKUS918IzC15jA0pdxH9OPBAyWNsQrnPxALg4ZLHKHtefwh4pMT9pwAblzzGA8TrVdR0YIOSx7iXSN9S1ErE56qX7qHcAJBViJRvZfyLGCVe1OrAWiWPcVvJ+69VO05Ry4jnUcY6xOtV1FPE+1HGepRbuPFJ4L6Sx9gImFri/o8B80oeY1NgUon7z69tRQ0Bm5WpENFp8WiJ+08lXqu8lxG/3WW9EvgSo8+n84jXt6hpbRx7LnF+Kmom8Tks427iPFvUqsR5vIyyvzuziOuRooaJ650yyv7uLKXxGmssc4CVS9y/nd+d9YEZJe7/BHB/yWOU/d15lPjOlrEZ8dvQK48Q1y9l9Pr6aDJxnVdGFddHZX932rk+ugtYUuL+7fzu3EG0U4pao7YV1c7vzmxgtRL3X0K8VmWU/d1ZRFwTl7EB8dnKWgl4W8lyIH6LXwNcnNnXzvVR2XZZO9dHZdtlRXTr9W+lnd/9jRk73lD290hqS3amhJubm5vb+NrKXLAAfL5k+WUbeADfK3mMf7ZxjLNKHuMPbRzjqpLH+GEbx7i/5DG+WLL8ySXLHwbeX/IYG7ZxjNeWPMaubRxjj5LHeFkbxyi77VyyTq9r4xjrlzzGB9s4RtmOpS+XLL9sABHgxyWP8fc2jnFeyWOc2cYxrit5jO+0cYz5JY/xnyXLn1my/GHg3SWPsUUbxyi7vaVknZ7exjFeUvIYL27jGDuUPMZhbRyjbGfix0uWX6Zzs+4bJY9RNiAB8LOSx7i0jWNcUPIYv2rjGDeXPMbxbRzjiZLHKLt9uY06lT3Gh0qWv14bx3h9yWPs3MYx9ip5jAPbOMY2JY/x9jaOUSboCPDpkuW30xl7fMlj3NzGMX5V8hgXtHGMv5Y8RtntrDbqdE3JY3y/jWM82KXn1+nrf2nJY/ysjWP8q0C5b4Ryo1ikdlxPdHLc1O+KSJIkSZIkSZL6y/RN6rWra9vTapskSZIkSZIkaQXVy7x7UlbZ3HSSpHKOAD7RweN/Qkw9zrqdmF5ZlGtKFOeaEsW4pkRxrilRjGtKFOeaEqO9DPhayWMfRKTNynNNieJcU6IY15ToHdeUKM41JYpp53dnIq8pMRO4nHLXQd8Bjknsd02J1lxTQpIkdd2WxMVzuzkpy+aeliRpRbIK0dFR9Lx6GQ4ClCSpiO9S/Py6hPJri0iSJKmHTqa9gIQdJ5IktTYEnEbxc+sX+lNNSZLGnYOIWXFFzq83UG4GqyRJknpsQ2IKZ5mAxBPATv2orCRJ48h7KR/0P7AvNZUkafzYiEi7VOb8+qO+1FSSJElN7UTkaCwakLDDRJKksc0E/kDxDpMTcBaiJElFHEnxmRLXU369E0mSJFVgE+BsWl/MXQ7s2K8KSpI0Dk0BPk/rEZ0PAv/WrwpKkjRO7Q/cSvPz6zLgJGD1PtVPXebIDUmSJq7nAK8BngWsRnSiXAv8EjiXuLiTJEnlrAccCuxV+/8y4G5iQMCpwKP9q5okSePWNODg2rYFMB14CLgEOIVYT0KSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJEmSJKmfJgOvBYb6XRFJWsFMB16Fv7+SJPXKBsCe/a6EJEkTxNOBZ/a7EpLGtyHgQODvwDDwnv5WR5JWOF8kfn/PB/boc10kSZpIVgf+A3gUeADYqL/VkSRp3JtJ9CEuAY7Hc6ukNs0AbiM6xIaBxcCb+lojSVpxfARYzshv8KXApL7WSJKkieMARs6xw8A/gU37WSFJksax1YDfMfrc+om+1kjSuLYPEeHM/qicCmzRz0pJ0gT2QuA8Rv/uPg7s0M9KSZI0Af2I0efbBcD7iZGekiRpbDOAtwL/YvQ59QoiHbEmIPNLqyqHA99h9Ajd5cA/gDOBK4GHianPUpUeBZYBdwJLOyhnBrA+MAVYtQv1ksqYBcwGdgUOBrbK3b6IWFfizIrrJUnSRDedOL++KLf/ceBPwB+Am4HHiPOxpO5bRrTrFgH3dljWbGK09krYGSr1yhTie7YJsSbT/sCaufvcTAy2m1tpzVQZgxKq0gHAicDa/a6IlLAYuIUIlP0G+DXwZIv7zwJeSXyutwc2J06s0qC5AzgU+Guf6yFJ0kQ1Dfg68C5sY0v9Np9o150PnA5cNMb9tyPadfsATwPW7WntJBVxFvBGYvCyJHXFbOA4YhTDsJvbAG8PAZ8iZkBkrQZ8hRjt1u86urm12uYBRwOrIEmSqrAXjekT3dzc+rtdBexHox2Acxi9/pqbm1t/t6uJYIRWAI7iUL+sRqQYOYhIN7IJLr6qwXQz8HLgWmLq4Ck4ekaDaSlwO3AxcEZtM02EJEnV2xY4hEhHsT2NKSkkVe+nwBHAU8DngKNwprvUb08ANwDnEhkrLu5vdVQlgxIaFFOJi/WV+10RrXDWIKbcbwE8g5F0TFkLiCn5n6Axr+jlxNTC64g0OYuBR3pXXSnpMWJq67J+V0SSJDVYmWjrTO13RaQJahrxPZtFpGB6DjEAMh8QvAi4FXhzbv9Cok13EXA98AAxuGdh76osrdCGiewUrisrSVLGC4BLaD2t8Exgx35VUJIkSZLU1EzgnUTHZ7M23SLgs8Ba/amiJEmSNNok4Ks0XrguBz7Yx3pJkiRJkopZD7iSxnbdXOBZfayXJEmS1NRPGH3x+oX+VkeSJEmSVMIcIghRb9MtAXbua40kSZKkFtYm1pQYJtaMmNnX2kiSJEmSyjqMkaDEN/pcF0mSJGlMnyIuXg/rd0UkSZIkSaUNAVcBj+AaEpIkSRoHVgEuBCb3uyKSJEmSpLa8Evh8vyshSZIkFbVGvysgSZIkSWrbJGD1fldCkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkiRJkqSShvpdAUmSJEmSpBXYbGCT3L7rgIW1/+8ITM7c9hBwR++rNdA2BdbK/L0MuKo/VZEkSZImhrWAzTPbbGClvtZIkiRJUi8cAQzntmdmbn8kd9uJVVdwAJ3I6Nfk0f5WR5JUxpR+V0ATzjbAWwrcbxEx6mN+7d8bgRvwQqKV2cBvEvv3AZ6suC7jxQbAkQXuNwwsyPz/HuA24nP5cG+qNspM4MXAS4DdgK1r+1JuB64G/g78GvhHBfXrtR8BW+T2fRL4Ux/qMl4cQLxGWdcA7+xDXSRJ0sS2K/CKHpV9M/D9HpUtSZIGlEEJdduWwMc6ePw1REfrScCt3ajQBDINeF5iv9/j5ubQ2edxOXAF8Cvgm8QIpW5aG/gg8C5g9YKP2ay2HQL8J9GQOxX4FnB3l+tXlR2Bp+f2ze5HRcaRdUj/HoxlDSKQmXU5Tv+XJEnNPQv4eI/KPheDEpImplnE4MMs216S1CMH0jjttJ3tKeArwCrVVn+grU/6tVqtn5UacDvRnc/jMDFj4r1dqtdk4KM0TsPuZFtMzDjYsUt1rNI1ND6fV/e1RoPvrTS+ZhcXeFzqO/HWHtVRkiRNDO+ge9es+e0PFT6PQWb6pvJM36RB9yxse0lNTep3BaQmpgEfIlLUbN3nukgQI8yPBX4CzOignLWAc4AvMnZAaQmR4my4QLlTgTcSMzskSZIkSZKkgWTaF1XhSmKKWtYqRCfqasQ6FBsDQ4nHbkXkld+VyPO/IquPhk/tV3E3Aucn9tc/k+sQC0tv1OTxryfW8Hh7G8feiJii/rQmt18L/BY4u1bP+4n3dyqRimpjYrTF3sSaAtPbqMMgWkzjZ3tZPyoyjiyj8TVL/T5IkiT1wi3AvV0o5+oulCFJkqQVXCp90+cLPG5dIq/+vxKPHyYCE5N7UF9NbKlUNd8r+NinAd8gFmVPfSZfW7IuqwBXNSnrllp5qcBcM7OIdFK35cpaWrJeWrGYvkmSJJWVSt/0nr7WaOIxfVN5pm/SoDN9k9SC6Zs0KO4jFhLeDjg9cfuewGFVVkgrvJuA9wMvJGYs5P0X5YII/8fohkXdb2v7f0a5WS8LgP8hvjOfxlHykiRJkiRJGgcMSmjQPAG8Crgscdu/VVwXCeCvwDsT+7ciHWRIORh4c2L/icDLic99uxYRAZLdiUCKJEmSJEmSNLBcU0KDaBnwAeCi3P5diHz/DxQoY3VgVWLB7OXEdNfH6E5qm9nEWhgLgQeJxYjbtTKx8PFQraxOOqcH2WxiTYQZwMPEc32srzUq59dEiqUtc/t3JlIytTIV+HJi/5VEyrJurZ1wBZGa5zNdKq8TU4l0VcuYONOo62uOPAk81ee6SJIkrejWItoXM4H5RPuiG9edqwPrAyvVyptHzFDuxHRgQyL96hNEe/bhDstsZRqRHnnt2vHuo7PnMIVoW69R+/upWrmdvi55M4n24lPENfeKZoj4jEC8tq4dWdwkoq9ondr/H6pt3fwcTQXWI75Xi4jfnFRGhSqsRvSxLK7VY1Gf6iF1xKCEBtUlxPoSm2T2DRGparJBiVnEIti71LatgM2IDsS8xcD1wD+A3wG/Jy5gW5kE7EWMdH8BsAONa1vcTXRM/x34c63uzU4KmwKHAnsQI9tXy93+eK1+VxFBmT8SF5F1hxMXmHV3Aj/OlbEe8L7cvjuB45vUqahDgR1z+75MnOzz1idmvLyUSH80I3Gf64HzgJOAv3VYtyr8mcagxAYFHvcyGhe2Xgy8ju53bj8BfLTLZRaxG7Hw9v7Ea7Rq7vYbgb8Qi3yfQXzOs/aslZH1JcYOIr6N0e/JA8DXmtx3G+DpxGLj04jX/gbgQhqDgbOAgxj5/M4mLkLrniBm0FwI/Ib47mdtSuPMrp8S3+0yVgU+mdv3e2KNnZWI9GJZ59MYzIX4/B1e+/+cxO2vArYeoy5zgWPHuI8kSVKvzCGuWQ4grh1nJu5zI9G++CFxrVbErFq5+wMvpvE6FqJdeh5wCvCHguVuCLyxVvbOidvvIq6Pf0UMgOrU5sDriVnYOyVuv4JoN36H1h21mxHt612Jtt/mxPXz1MR95wPXAJcT1/h/ofggwGnAfoy87usz+j1dSlxjXwScRbSLuzWYK2t14CM09o39k8Z2drdtQ/Qz7E/MwF8jd/u9RHvjfOJzcl/u9oNpbEMdRflgxm61srK+QnR270m8T1lHkx5guA7wbKIttDIxOPR24vOfTXt8LfCjxOOHiM/cNkQ7ezLRbryBaIvn3//tgFcSbbZdaeynWUb0M5wLfJv4HmetT6zRmPVz4ruSvU/9e/VcGrPNXAOcSqRVrgdFe9H22oV4ri+o/X967v4PEn1IVwIXEJ+Z8TQIVJK6ot2FrlPOSpSVXVz4RcSJLrVwcJHtMeKknf9Br3s1cRItW+6TRFqebAf+ZkQgZFkb5V0CvIG4ULoyd9uFiXoPERfk2fstIU6o7ZpBBB+yZeY7YiEupI4hOm2LPr/lxEVWkQ7+sjpZ6Drv6ERZRT7bv0887qQ26zBoXkB8Bst+P04Gnpcp55jE/VYqcPw/5B5zbe72VYBPEQGwZvV5Tub+KxNBgPkln9NNRDCo3pB4fuI+ZRdGhwhA5supBynmJG77VJNy9iv5fFLblW3UX5IkTQz9XOh6dSJV6eOJOrTafsvoAW55M4EP09jGGWu7AHjGGPU9hphVX7TMf9Xqm99fZKHrOcTadYsLHusWYtBdyv+VfC3y2w00dmDnTQLeBNxWsuz64KPNM2WdmLtP2Zkya9DYvh4mUknPavG4Tj2NCHCV6ctYSrR9Xs5Ix/ixifuVWfOw7shEOVvVbjsqcdt6mcdOIvoqLqjVMVX3Bbm/f5k7/hzgq0Sgrtnzz74fmxFBjTJ9K08BxzG6jdlqEeo1gS8Sbdci5T8I7F17bDfbXs8ngg1lH7+ECHbW6yRJK4RuBiV+kigrm5e/Gz/2w8TIjtmZcicD3+1CuefUyjuIuEDqtLwNKBaUgBjtkX/8UU3uW8ShifLelbvPdsAdHTy/u0mP6ulEN4MS30iU9eExHjObxouz5cD2bdZhUEwGvk7nn+k31crrdlBiiJipcH+BOtSDEtsCt3b4fObWyjIoIUmSJpJ+BSW2IjrR271+uY/RA1Dq1iVGULdb7qPE6Pa8zen8ejK7jRWUeG1if5FtKelr0293qd7HJMqGCAKc12HZS4EtauWdmLutTFBiVWLwX778XgckDic6yDt5Df6nVla/gxJ7UqzDvFlQYhrwWWKw6Fhl1N+Tlxa8f7Ptskz9mwUl9iHadWXLXkjMsu9W2+somgd6im53IQ0w0zdpkK2e2DevxOPruSjruSmb2QU4jYgiDxOj2N/Y5L7zattC4sS4DukpvnUvBX5BerprNp/oTKIDezbdWYD+RGJEUXYWyNsY6fwt6/Dc348xesrlbsSU3fyUU4iRP39jZBTSDKLTdvPc/TYgZhU8k5iqOmjWTewb6yS/D43TSP9B44j+8WQGMa31gBb3eYR4v9eldXCh1feyXesSAc0yo0L2ImbrpH5zIEbh1KdLr0fz72gvno8kSdKKaBdi5vzsxG13ESlWHiIG/EwnUuPmU63OAc4mZrDXU7dsSQxs2TRR7v1E2qcHic7AKUQ6mqfn7rcqce24GyOzx9evlZtv4/TKS4HDaK8TejLR5r2Fcml0HydGYK9GYxsn62PEgLP/zezblGjr5dPaZt1PtOHXI91+pnbcZrcVtQrx2Xpubv/lwL50f62MumOI16aZhUSbYxbpdnVds0wPVZlKBESOpL3PH8T38DQa00O3chgxeDTfjzlMtK+vJdqhEK/hXsT6D1m7ErOS8mmv6t5FpJtqp09mBpE67gNtPDZvDvDfTfgbMA8AACAASURBVG57hPicPE58lteh9edFGlgGJTTItk3sm5vYB/GjfB6Ra/0aIn3RI5nbZxI/7M8jRvrsmXv8nsRJ7mEaAxJPElNZTyTSwORtT+TBPISIjNdNBk6g8aLpbGJ9h7NozLk5m+hM3Q94DZFOph0PEhfKh2b2bV6r359LlrUJkSor6xRG1gWYQwRe8ifCnxEpj65pUu7ziHUDnp/ZN5sIduxDe8GTXnp+Yt9lYzwmf6ELkd9xPPtf0gGJS4jvyTmMDh6uS6yh8hLiM51a76Vb1gYuJqb0FrUBcXGbD0gsBL5JfLb/ykgO02mMfOfruUUH3VXEaw/x2nwxd/vxjP270KvGmSRJUt5soi2TD0j8kpiF32wG567EdU42PdEsYsDKHkRn7i9oDEj8kRix/RfSbZBtiXbNyzP7phMpSZ9FrCf4Y9IBiQuImdrnMDrYsVXtsXsRueLLdiquk/t7IdG+PINoC99HXIfvTLRvn527/wzi2n03IrCTNUysqXAu8ZpcTwQwlmTuszrxfF9EdE5vnCvjv4n3616iLX46jQGJYaLN+CPiWrS+1sWkTNkvJbIxdGPgHpm67J7b3+uAxDtJByRuI2ag/672/7rVifbyvsSslk5SMXfbaTR+nsp4BtG+WrPEY3Yh1kPJ9mEuJtqmXwPuSTxmErG+45cZmV0D8Z04mljzL2/X3N8LiDUEzwJuJtL8rkW0Ad9GrDmatTHxve6k7bUh0U+StYwIJH4XuDTxmE2J9ukBxPelVdBQkiasbqVv2i5RznxGn4R2JU5Cz6fcRcoQcdGZL/9S4uSY3fcojSeaVp4J/IA4aaTy2H+2RFmziPRA9RQ0ZdI3QVzg5o+fWlBqLJ9JlFMf0TBEYwqdxYyk5RnLNGLUfb78fduoZ0q30jftmyinSDqb1PTkQ1s+YrC9hcbns5S4yC4ySmY2o3MCv722v1vpm/Lbk8Tn6w3EZ7Y+G2lz4vvxXiJoln/c1TQ2rJrZkwg01X+jYPDSN2WlvhNvbfkISZK0okulb3of0THf7tbq2vGM3LGWEh2ARUwhnQb4YGLASX7/f4xRl7ohYvHf/OPfQ3Q25/cvpPh1/0wiJU++jFbpm+rbfcDHidkLrer+QdJrGLw+c7/31sratGC962aRXmfuI7Xbv5e47UEaAwPNbEOMQK/Xf5va/hNzZY6VvmkG6fX+ep2yaWcicJU/7nEUm/UwA3g30fE+TAx8hP6lb8puy4j3/r1EEGU9IsC2KdFf81biM9qqjIeJDvdXEIO/1iSCblsRA9uOojFN9AM0n+2QN5to32Uf/xTpNn59e4joi5mZKK9uMunP9q2Z+7Rat6KZk2h8jQ8s8DzrNiXW6FiI6ZskrWC6FZQ4JVHOyV2qY12+QzO1UNJX2ix7Zxrzn95Le7OTZhHR/dmUC0qkFrxeSLlROJNovAC4PHN76v0+okT5ENOf78yV8fuSZTTTjaDEasRIoXw57yjw2FQO3PG6nsQM0rk12+nQ3hK4iN4FJRbWyhzrs/72xGNvJEa/lDFENBTurv1tUEKSJE0kqaBEp9vWTY61T+K+/16yvivReB1+OY3tvfwI5rFMpnEdgtuITtV8u7JVqtOUI2h83mMFJU6kdTAi73OJMv5Ysp7NrEvj4t7nEwP88sGQxxn93Irah+hkbScoMZWYIZF//n+ltwEJSLdZjm+jnDWIgNugBCVOpfn3OCvf1q9v84ng5rQxHv/p3OMWUn62+jY0Lgb//Sb1+inF+0ymE23AfBn1AW5lgxJDNH7XzyhYl7zNiOCEJK0wuhGUeFuijGGKR8KLemOT42S3tzd99Njyi5z9pZPK1pQJSgB8iMbn9O4Sx0s1CrKjlM7P3XZWibKz8gtzL6E7aX46DUqsT0y5zpdxBcWmRD6ZeOwgTb0t492kL0TbNZmRzv9uBiVuI6b3Fjn+zbnHLiVG+LRrTu1fgxKSJGkiqTIokR/J/uc265zqaM1uN9De2gSHjFHuMJFOpqyyQYmftXGMGYyMtq9vyyg+Q3gs+U7/+4mBhfnndWQHx1idkXXcTsyV2ywoMYVI25WvRxUBiV0Tx72eztaiq7c5+hWUeIKYQV9UKihxGcVS7s4kZkVkH9vOoFeI70y2nNSi9P/VRrmptuwrareVDUqsnbh/O78n0rjQrbx8UjdMIdIbnZC47UwiV3w33VjgPmVy0+fl14PopKx2nURMFc0qM5Mhv8D1AkZyL25L5GbNarYY01hOz/1dX1SuX7YgRmT8A3hB7ra5wKsYWWOgmUmkp3vOT+wbD/IBumwHeTuWEdNiu+leotP9igL33YvGxRB/Q4x+a9f9HTxWkiRpRbc5kRc9q1vti7wvM3qNhKLOIlK/NLOI9jtNy3hy7Ls0WESkQcqaRON6i+26Iff3bODVuX13094sgbpHaGzftjKZaL++Irf/MiI1UK/XTUsNcvwc5Z5DXr/bHC8i+hnadSkxiOv2Avd9JaMXq36K9rNZ5H8TNk3c57bEvrGk1plcr41yIL2maD/6kaRKGJRQv61D5PL7PJEm6DM0fi7vongO0SLWIDoji4wIeRXtj9i/L/f3+nRvrYSiHiIWicvaiYjYj2VNRi/mBjFdtH4BvFfutvuJlDztuJnRC5NDsTq2Yz9ihH1+u5AIQjxOTPf+HI1pfG4hnvcdBY6Tyg+6hJhuOt6sQePaKucT79sgWUDxhsULE/u+08W6SJIkqZw9GT3Sez7tpxe6ixhhnbKcWIS5HU8RqV2b+TWRzmlQnZ3Yt2NiX1FTiZH729KYSmoSjbPLT2TswV3dMolYGDgfGLmM3i5qnZVvc6Ta5+PNvzp8/FwilVIR+T6HP9L++5YfuNat/tB7E/vKLlxfN4/G78eLgI3aLE8aaO3kt5fK+hCNKYMmEVMvxzIXOIjGDv5mhogLomcSF1fbEx3La2a2ViefhYwe3b4VcdHycSKX3/KC9YDooM/nyvwVEdn/OtVcBEHMPHldbt8RjJ2b9fU0TivNzmLJX2D9nXKvT948Rn8myub1L2p92kuh9EPgAxRvZKQutKYSwYpWo6sGUWox+T/1oyJdlJ/lU1+kTZIkSWM7luiAb9ediX359sVVRHrNdj1ADILLu43OZi83C3YAnNdBuVW4KrGv1foOM4nBYvX29WZEO63exl615PGrakMMAd+mMc1QPSCRHxDXC+sBT8vtu5jx1xbspz1zfxeZEd9Mq+9tJ1KfpXbTcz1BfEd3zuxblZhd8mngx/j50QRiUEJVmEF7P8p/Al4DPFjgvjsAb6jdf9M2jlX3FxqnDG9LXHDPJXJRXkycKG6i9SiP7wD/xujO3JWIk8nHgHOInKlXESP0H+ug3q2cD1wHbJfZ90bgo7Se9ptP3XQxcE3m7/yCzS+is1FB+ZE17Y4u6LYzgKOJC4EylhGBifzCXWtQPMg2KFI5f/9WeS26K/+cbiRmyUiSJGlsN9P9DuZ8+2J3OmtfNOsw37TDclvNpL+yg3Kr8Agxw31OZl8+cDOZWKj7tcDBdGetv7q/d7GsVr5BY9riKgMS0BiQgPHfhqrSTCKlW9aHKL/wfa91e7Dpt2lMab4+MevnK0Tf1PlEP9K1FJ91Ig0cgxIaRHOJxYJOYOwo8CbEBcfLaG9Rp7wziIuHTRK3rUecAOsnwSeJi86/ESOsz2Z0p+ZVwHHA+xJlTSdmgBxU+3s5cWF/BXGxdAbt5TNs5nvAVzN/r07k1fxxk/s/i8b0SfkT45q5v6fR2AHfiSILHVdhGu2PyHgA2DC3by3GX1AiFSCaW3ktumeIxpk4qWm3kiRJqk6v2xd1U+jdAKi7elRuN81jdFAiOzhsf2JWf7OFyDvxBNUEBD4DvDe372oi0FJVQALSM//HcxuqavnfA4g+gkHpJ6h7osvlfQ94MxGUzZsFHFbbINJDX0P0V1wE/I5ig3qlgWBQQoNgObHI0d+IhcN+RrG8+3sRixXlR9jnLSZO/nOJH+j6Ngn4YO6+TwJ7Ewve5kfq5K1EnCh2Jy56Ftbq/gViFgVEJH8Kkb6qVdBkEnHhtzUx4+NYYsHdrwM/JxYV7sQPiNH+2RkrR9A8KJGfJfEwcFpuX6/SK/Xa2YwO0GQ9n7iIzdoX+BqNF7ZF3EVjUGIjYkTDeJK6IKzygr7bVqXx/PdoPyoiSZKk/2+8ti+yxsM1Zb4Ttd6ePhr4RIHHP060c+rt6geIQMezab2GYhWvzceBzyb230j1a31MtDZU1cbL70GnfTV5y4ADiSwde49x36nEmqE7Ae+oPfYM4ItEf5I00AxKqArHk+4EXkRcmLSTMqWeUik1JfcOIv//JcRshWaj0rehMSgBMUNhJ+CdxKyI1LTLlJlEzsrXE2sPHE+cFP6dCFZ8BHgpjYt9NfO82nYWMXW2k/RODxEntTdk9r2QWPD7ltx9ZxLpnbJ+yOhA0SQafz+uJU6A3dKrqb33Auc2ue1cYhbJ+3P7jwTuBr5U8lg3EO9h1s6kF5gbZFMT+6paoK4XJtrzkSRJmgjysyJuoLN1K+reweiZEfcCP+pCua8CtsjtG4+pVBYTs/tTAYlhIlXMz4nR2NfSvP3+KVoHJTpZH6SIVYgBgimvIdrE+bUueynV32abo7jULKlz6V4KrDk0rjkyKBYA+wCHEn0Tzy74uMlEFpGDiYGVH6WzdT+lnjIooSrMp7upiAC+RWNAYj7ReXwynf/wLiZSLx1HpDE6kJgRsVviuHnTgP+rlfHd2r4La9vaxAnihcSo/M0K1GV/4mK81QVeEScwOigxBLwV+GTufocQ0wLrhmlM3bScOFFmR3/cRoxMGe8+RDQuDsrtP4YYEXRyibIuI17jrAOJUUjjSWpUUZGF6gfVfOIznF3vZTw/H0mSpIngYSJ3et1ddKd9sSfwnMzfj3ap3I1pDEqsRO/WCuyWWbm/HyHdmX8FMYP+msRt3Thut2UzEywmBiJl972L+Ix9qsf1qJtobaiqpV6/C4H/7FL5z2JwgxIQ/TAn17YtiGDDHkSf1NpjPHaI6Neo/ysNpElj30UaODsQP8ZZDwC7AD+h+5HgK4H/AvYjRtg8E3g7EXC4tcXjjqVx0bB5jOQI3BzYAHg58N/ECbbZ6JG9iBE+nbiQxrRBh9MYnMynbrqAGKWUl79ImN1+1QbKciJ484/c/iHivcvPfGjl/MS+Z1N89s2gSF0Qjuf3ezmNU6d7/Xym97h8SZKk8a5X7YtelZtKRzQeOp3zaYWmE7Pls84gOj+7FZCAmMkwY8x7de5hYgDgkYnbPkmxFFXdqkeebY7iJlobtBO3EjMfDiH6mLYG3gT8D9Fv0SyF1AeBXauooNQOgxIaj/ZJ7Psc3Z+NkbKM+NH/LhGY2BLYCvg+jcGQlWhMg5R3L7EuxieJQMuawL8RwYu8ToMSMDJzo25dYiZG3ebAi3L3yc+SqMsvDPx00mlxxqPHiBkt9+f2zyTWG9myYDnXEwuYZ00iUnmNJ6kUaM+svBbdlX9O29Hbz+9Ya9/000T53kqSpPEt377YlsbO8m6UO5uY5dCp1IKyRdsJ/TKHxlHW+Y7spUSbdEmXjz0EPKPLZebNI/oLLieyF6RmgBxNpKvqtX60oQa5zVHWIzSuf7JTPyrSZd1oe91ErA/6PuIztT4x0PWpxH270Y8k9YRBCY1H2yT2/a7yWoy4hVg0Op8GCWJ0SRmPEUGA3Wg8Ae9A59/Zk4jFvLOOyPz/zYye4voAsRZFyoW5v1el3CyCQfcv4BXE2idZs4nARNHpxz9M7HsT3WkIVeWixL7x/l7/Jff3DGIKbzekGnBjpX3rp4k0okqSJI1fqeuzF/SgXIAXd6Hc1CyCXne6d2qXxL5Vcn9fDdzTo+P3sg2xhAhIXJnZ90nS64d8hRh13kvXEWljs57L6PZ2J1JtjokUlBim8bv7HMbHbKRWetH2uo/4rB+auK1sn5RUGYMSGo/WSOxLjVKp2rE0Lly1bptl3QL8NrdviM6j6gtoDDIcAKxHLIp0RO62H9B8sbbzEvve21HtBs/FwGE0TofcFvgV6cW38r5LY2BjOvBTuj9CfYjeLN52M7GAfNbBjJ3LcpClPr/59T/alU8NBY3T5PslH+wEgxKSJGkwnJvY1432Rarc99F5f8jfE/vy69INmtcl9uVno/SybZ1vb3bTIhpT8A4T1/in5fZPAU4FXtrD+iwH/pjbtxmxxkk3pNKHDUqbo1vybbYpVLtYeadSba9epjA7ncYMIu32SUk9Z1BC41F+pD8MRufoQmJmQVbqJFTUHYl93Vgv49u5v6cQMyT2BjbMHatZ6iaAS2hMb/QKYvTHRHIq6YWp9wS+WeDx9xG5HvN2B75O90bKrEXM4EgdqxvOyv09jc4bifn1TKp0Ho3BotcT67y0q/58HkrctnMb5fUirVKq8eKFqiRJGgRXAHfn9h0AvLDDcu+ulZ31DNId9GXcCszN7XshsEmH5fbKejTODriUaMdmdTNv//W5v3cA9u2gvEnEYLoylhFplc/J7Z8K/JzuzMZpJpXRodNUvq3aHO2kNxrkVK6/oXHg54cZjP6fIsq2vbrRPr4j93cnfVJSTxmU0HiUzwkK0aFexsbAcYn9hxGdu+1YjcYLuHynZxmb5/5eSHfyev4F+Gdu3+E0jlr5I60X8l5ETHvNGiJGoazfQf02Ji4Y1+ugjG77NHByYv/hwMcLPP5oIh1U3nuIlFpFZly0sh8xUuvADstp5Ws0XhB+jPY626cD3wDe0mGdOvEgcGJu32pE0K5sQwfidainuZoH3Jm7fQ9i/Zmi9qExRVo3zKOx4Vn291OSJKkXlgBfTOw/hc5Sn25AeiDQt+gsx//qNC7GO43I7V7G7sR1dS8NEde5K+f2f53G9vUOlA9M7EesQ5F3bGLfCaSzH4xlPeD3lLumrlsMvAr4W27/TCJDQbfSuOb9hMbXd3+iHVlWfVb8V2t/558LJcudBhxT2wbVjTTOclkT+CWdzfZ+NsUGGHbqQRrbXnu1uP9ewJk09gUVNUTMxsm6o82yJGncOZCYIpndPt/lY+ybOMYdxKJdY5lORNYfSZRR3+YBH6Dxgm0sn02UtZxYU2DrkmVtS5y8smX9jMiPmd3XbqflBxJ1XZr7+7UFylmZuMjKl3Urke+xjPWBLxPragzTnaDETom6fa/NsmYSI4lS7/HrCzx+D2LhqdRn7lLg+W3U6fnEdPT8+9grP6Wx7g9RLk/lIUQ6qGFisXiIC+F8uSsVKOsPucdcW6IeEBdsixPHPpXi02rXIgKcSxidM/aURLlnM3bAY3PiN2NZ4vHDjKxdMydx26cK1vmSxGP3b/kISZK0InsHjdcO7+nRsWYQg3lS7b3dS5a1LnGdWW/7/SNR7jzKp1yaBXyCmCU/TPq67aOMPSN6IyJQ0Oy6Lxswybdff0nxQaaTgO8kyr+JGJl9QuK2Hxaof/05/IhoE6Wew2winUx+/1UUb+9NAz5EXGsPM7LG5Im5MlOj0vNmE7M38vW5n/Jt9qI+mDjeEsoN0NqFaPsPM5LNYAYxCj5fdpGFjV9KDGpLvWfDjAR+jkrcVradfmfu8b8s+fjtSLfZLqR8sPIZxOd1CSPf3+zWTjrfVRPlZPvALk7c3qzt9eLa7U8SAxvLzgh5beJYHypZhiSNW1UEJaYRF6X549xFdHKmFj56BjHaPXVB1Gx7EPhPYFdadySuTaTMaXYhVr9QPY3Iw99qYaZJwBuIi6L8459L94ISs0hfwNS3eyg+dfDFNAY06hdaJxPTYZuNYliPGM1xGjHzopOLnZRuBiUgGjb5i6phIoBUZNG215B+rerbWcQIo2YXxLOIoNwxxGyXVBm9DEpsQAQhUu/1scCWTR63LvA24PLc4/odlID4XUi9jrcAryYdnJxEvN9fZ3QDMRuUeDHpcv9EpP7KfifWJ9IH/IL0BXd260ZQ4j8Sj51PLIyW/66uTPxuFZkRJEmSJqYqgxIQKZCWJI65lGg37EnzASRziNnvp9A4yGtfRgZAZbflxOCR/Wlc9LluDaLD7yQaAwS/SZQ5DPyZuI7KpsdZjWgz/zBRv/zWKigxTLRtP0njyOi6mUR63asSj13CSPtl9ybHP5eYUZtP7zMNeAkRUHlyjOdQb7+k2srzibU9mnW8bkN0iueDVJ0EJSBSFqcCX3fS2YycZqbRPABwOhFwSAWYViXevzMYHbjKplhOBZuWAl8iOvPrphIDHz9CzLBo9Z4NM1hBCWr1TtXzMWIG/g4077PZmhiUeS6jP4dVBSU+lbh9PtH+y7e9Dsrd7wliANyeiftmrUy8V/k+lXl0Nx2bJA20KoISECPTW51E7yZG699D806+sTr/stsC4iT2A2JK8XHEyfQK0hfMrS4wl9YedxpxAv0KEa0/j5HRH/ntS7Xn3a2gBMQFdbM6ln3P3kHzET7DxMnxSuKC6hfESO17aR3IGcSgBMCOwOOJcufRvFM+6xDGboAME5+rucTn+FZaz+7Jf756aU9aNz5uJz7LpxJpuFIX/PVtEIISk4hZSK0+u5cDvyYWN7+U9Ptfv7isGyI9IyG7PUrr4GDq96AbQYl1Wxz3cWL02BVE4Kv++3ZlwbIlSdLEU3VQAiKw0Gowz2LgaiLVyc+JNJr30Lp9sRmRZqjVtexS4prybOJ69kKiU7VVXQ4k1pFrdd23gNbXfam26VhBifq2nGj/Xkhcr/6SuA5t9Tw/kXu9f9Hivk8RAZBbiYF7ZZ7DrFr5H27xuGXEdeeZRBv5fNIDoepbp0EJgO2bHOOftJ/OuZVNiEGUzZ7T/cAFxPM/g7geb/aZywYltqZ123IZkWIs1WdR31JtjkELSgAc3+I51N//vxLtttOJ4Eurz2tVQYk5FG97tfqdeZJIxf1jog/pG0Tw9SLSn4HlxKBISVphVBWUAPi/xLGKbn8kRne0+/ixToYHdLG8/2Fk6mw3gxK7NTneUtobIfIC0qmc2t0GNSgB8f6mLhiup1h+1O0pNkKl7PYnGhev64WdaR1sKLoNQlACYlTNMbRuxBbZ5ufK3ZLiwaT8djoxwyu/vxtBCYgZSmXqY1BCkqQVVz+CEhAjyG9PHLvdrT6jYGvgmi6Wuwdxfd+qA7TV9j3i9czvLxqUKLMtI/3erQHc0GaZS4EvEFkG8rfNyhzjjYw9s6LI1o2gBMS6AqmZM3+l+YyZTqxPOpVP2e2EXLmpz06R7SliBsJ7E7cNYlAC4M20Du6V2aoKSlArtxt1LvOdfEsbz0WSxrUqgxKTiVkLZToSHyBOCEPExUz+9tNoPYJhrO0PwBa18vcnRoy329E5l8bIdjeDEpDOq3pmB+VtSEzjzU8bLLotJt6DF1M8R2orvQpKQCxGl3oOf6bYwtWTiY7hVE7TMts8RqasVmlt4oK4nff6SuBNjLxO/Q5K1B1I43esyLaQmDq9XWOR7EiMXCta1t3ExTbEjIb87d0KSkCkCis6Y8yghCRJK65+BSUgOkCPp/3O7CVEJ+h+jG5fzCLakgvaLHcZMSv4VYykN3o2sTBv0TLuI17bIeCIxO2tghKXAteVrPNcIj1pM+10ml9CtLkgnaYmG5Sg9pzOpHwbeTkxg2DvTFkn5u5TNigBsA/p9sx5dLaQcjPTiPSx7QSw7iSCCKl00EfSepR9fvsTI+3HIxO3D2pQAuLz9mtaZ2potT1CLHL9usRtvQpKQLm2V5HMCs22q4nfImngFc0ZLxV1L9GpnPXPHh1rGdExfAaRI31v0nkElxIzI35eq9uC2v7HaKzrt4hcoc8FXk4sJPxMIh9nyjCxvsV5RDqk+kUcxPoAZxGLfx1C5NLchehobGYBkSrmJ8QU2sdzt/+ZCJrUddoJ+3niQjrrWx2UdzdxYf1pIjK/B5GjdLUm9x8mRj9dBvwW+B0j7083zKfxPb6iS2V/ibggTKVs2pf4XLayDPh+bXse8T68hBi51eq3eSkxy+J84vNwHnFxU7V5wDuBzxCNqH2I55G6eF9GfFZ/T6wzku/gvo54/7OKpKG6gNGzE+4u8JhWzqhtLyHej72JhadT5hPf918Q0+SbfW6vIkbOfZB4veY0ud81xFTc/yUa3RANpPzn9/rMbfnX7KYmZTfzLeLzcxQRRE39Nj1IfD/PKlm2JEmaOO6g8brj9oqOPRd4N3HN+RZivYndaezszrqdaFOdQXSAP5y4zwKiLXk0MSDkRbVyWy0sew9xHX4m8XrMzd1+GdF2fD9xfdwsteu1RHvvOEbae7fTeN2Xvb48ndFt0u8RA3QOINZu24/mA6PuqB3vy0SHbDP3EqlaDycWx21W/4dr9fk5keaq3v69LvEcluT+vrpW5x2I131v4jVLDUhbSLzevyZSad2Zu/1yRq//trBJfVs5l3i+qdnmB9P4fDq1mJhR8lVikNb+xGe62bqTtxIBhJ/V/l3W5H7H1W7/DDHQKrXmymLi2v+bjP4+30zj86x/Lq9P3LaoSR2aOZPRKbEuK/n4vL8DLyPWyHgTkbFhV5oHkZYSs4AuJj5LfySew2Y0Prd2ftfqa91kpfrA6m2vjxPveyq4k217nU/0Se1FBGKa9alAfHf/QqxV83t6n9JZ6oqhse8ijRtrMXIhOZm4OLqZOGm1M2qibgrwtFr5axIXe48SeR9va6PsDYgZBWsS00IXESefe2i80JoIJhOBmVm1bQoxVfRx4vVrdWG8IppJfEbWrG2rEhfz84mLjbtp74K7CjMZ+WyvSlw0P0J8Dx/rY706sR6xONiaxCi4J4j8szcz0gArahIRmNycuKhcTAR3bqZ8QKEX6u/dVOJ37SHSjXhJkqR+msTo9sVU4rqq3r7oZJBTttzpROfeY0Q7bV7JsnYi0uKuTVzPP0hc993YQf2amU7M0N2YSMW0nKjvXYzMBC5rGyJgUO8wf4oYcHMt3e30XIMYHLMm0Z54imhDXE9jYzo3XwAAAO1JREFUUGMimky8b2sxEmybT7x3D7RR3srE4MC1ieDE47VyrqC7AwAHyQxips8axGs4RAz0epxoZ5UNpFRhA+I9L9L2GiKChOsQz3Em8dweIgJXD/W6spIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSZIkSYPq/wHSSUoWhggnAQAAAABJRU5ErkJggg=="
    }
   },
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Interconnecting elements\n",
    "\n",
    "Fig. 1 shows a general representation of a neuron as a parallel interconnection of a passive membrane, consisting of the membrane capacitance and passive (leak) resistance, and a multitude of *feedback currents*, i.e. elements emulating the action of ionic currents.\n",
    "\n",
    "<div>\n",
    "<img src=\"attachment:neuromorphic_circuit.png\" width=\"500\"/>\n",
    "<caption><center> <b>Fig. 1: </b>Neural circuit as a parallel interconnection of elements </center></caption>\n",
    "</div>\n",
    "\n",
    "At the moment, our neuron consists solely of a membrane capacitor, so let's start forming the neuron by interconnecting elements in parallel one by one. We first add a single passive resistor with conductance $g_l = 1$ in series with a battery of voltage $E_l = 0$, mirroring the *leak* conductance element of the classical Hodgkin-Huxley model $I_{l} = g_l (V - E_l)$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "R = neuron.add_conductance(1, E_rev = 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The method `add_conductance` takes two arguments: the first defining the maximal conductance, and the second (optinal) argument deifning the reversal potential. Ignoring the second argument sets $E_{rev} = 0$ so we didn't have to explicitly write it here."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can check that neuron consists of a single resistor interconnected with the membrane capacitor by measuring its current-voltage (**IV**) curve. The method `IV_ss` sums up the steady-state IV curves of all elements in the circuit (note that at steady state the current going through the capacitor is $0$, so we can ignore it). Since so far we only have a single resistor element, the function should measure its IV characteristic."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = np.arange(-3, 3, 0.1)\n",
    "I = neuron.IV_ss(V)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Neuron steady-state IV curve\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V, I)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The neuron's IV curve is a straight line - characteristic of a linear resistor!\n",
    "\n",
    "The neuron is now a simple leaky integrator that integrates the applied current. Let us check its behavior by applying a step input current at $t = 5$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trange = (0, 15) # time length of simulation\n",
    "i_app = lambda t: (t<5)*0 + (t>=5)*1 # Iapp(t)\n",
    "\n",
    "sol = neuron.simulate(trange, i_app) # call neuron's simulate functionality\n",
    "\n",
    "# Plot simulation\n",
    "plt.figure()\n",
    "plt.title(\"Passive membrane step response\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Membrane voltage\")\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Input current\")\n",
    "plt.plot(sol.t, i_app(sol.t))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The neuron acts as a simple first-order filter: the voltage follows a step change in the applied current with a time lag given by the membrane capacitance and the conductance of the resistive element.\n",
    "\n",
    "In order to build interesting dynamical behaviors, we need to interconnect *active* elements that emulate the dynamical properties of ionic currents. All active elements (currents $I_j$ in Fig. 1) will have the same simplified form:\n",
    "\n",
    "$$I_j = \\alpha_j \\tanh (V_j - V_{off,j})$$\n",
    "$$\\tau_j \\dot{V}_j = V - V_j$$\n",
    "\n",
    "The first-order filter defines the timescale of the element, so that increasing the parameter $\\tau_j$ renders the element slower. The other two parameters are $\\alpha_j$, defining the gain of the element, and $V_{off,j}$, defining the mid-point of the hyperbolic tangent function.\n",
    "\n",
    "Importantly, we can also change the **sign** of $\\alpha_j$:\n",
    "- $\\alpha_j > 0$ element has a purely **positive** resistance current-voltage characteristic\n",
    "- $\\alpha_j < 0$ element has a purely **negative** resistance current-voltage characteristic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Designing a bistable neuron\n",
    "\n",
    "Let's start by adding a single active element. Since with this element we would like to capture the fast excitatory nature of the activation of the sodium current, we can set it to be *instantaneous*, i.e. we can set $\\tau_x = 0$. We would also like this element to have a *negative resistance* characteristic, so let's start by setting $\\alpha = -1$. We will keep $V_{off}=0$ for now."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = -1\n",
    "voff = 0\n",
    "tau = 0\n",
    "\n",
    "i1 = neuron.add_current(a, voff, tau)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's first measure the element's current-voltage (IV) characteristic. We can do this by calling the element's `out` method that takes the input voltage and returns the output current:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = np.arange(-3, 3, 0.1)\n",
    "I = i1.out(V)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Negative conductance characteristic\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V,I)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interconnection of the element should have changed the neuron's steady-state IV curve. Let's measure this now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "I = neuron.IV_ss(V)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Neuron steady-state IV curve\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V, I)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We observe that the steady-state characteristic has indeed changed - the addition of the negative resistor has decreased the slope around the origin (note that the slope is exactly $0$ at $V = 0$).\n",
    "\n",
    "Let's now try varying the gain of the negative resistor in order to change IV characteristic: we wish to make the neuron **bistable** so that for some $I_{app}$ there are **two** stable equilibrium voltages. Let's change the gain of the element $\\alpha$ so that for a range of voltages there is a *negative slope* in the IV characteristic, i.e. the IV characteristic should resemble a letter \"N\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change the value a = ..\n",
    "a = -2\n",
    "\n",
    "i1.update_a(a)\n",
    "\n",
    "I = neuron.IV_ss(V)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Neuron steady-state IV curve\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V, I)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the IV characteristic has the \"N-shape\" the system will be bistable and we have created a *switch*: the circuit will remember its state even after the transient input is gone. Let us test this behavior by injecting a series of current pulses:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trange = (0, 100) # time length of simulation\n",
    "\n",
    "def pulse(t, t_start, t_length, mag):\n",
    "    I = (t>=t_start)*mag - (t>=(t_start+t_length))*mag\n",
    "    return I\n",
    "\n",
    "def i_app(t):\n",
    "    p1 = pulse(t, 10, 1, 3)\n",
    "    p2 = pulse(t, 30, 1, -3)\n",
    "    p3 = pulse(t, 60, 1, 3)\n",
    "    \n",
    "    # More pulses can be added here\n",
    "    pulses = [p1, p2, p3]\n",
    "    I = 0\n",
    "    for q in pulses:\n",
    "        I = I + q\n",
    "    return I\n",
    "\n",
    "sol = neuron.simulate(trange, i_app, max_step=1) # need to set max_step so that the pulse is not *missed*\n",
    "\n",
    "# Plot simulation\n",
    "plt.figure()\n",
    "plt.title(\"Switching behavior\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Membrane voltage\")\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Input current\")\n",
    "plt.plot(sol.t, i_app(sol.t))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can play around with the pulse timings, magnitudes as well as lengths to see how they affect the switching behavior. If the negative resistance region is increased, the system has a larger energy barrier between the two stable states, so we may notice that we need to increase the pulse duration and/or magnitude in order to reliably switch the state."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## From bistable behavior to spiking\n",
    "\n",
    "Having designed a neuron with two stable states we are now missing one more component to make the circuit spiking! This component needs to emulate the slow refractory nature of the activation of the potassium current. From the dynamical systems point of view, this current essentially needs to \"automatically\" adjust the applied current, so that when the voltage jumps to a \"high\" state, the slow current *decreases* the applied current and drives the neuron back to the rest.\n",
    "\n",
    "We now introduce a slow element to our neuron in a similar way as before. Unlike the first element, this one needs to provide *positive resistance* and needs to be non-instantaneous. In order to make the element *slow*, a good rule of thumb is that the timescale needs to be at least $20$ times the membrane time-constant. Since the membrane time-constant is given by $\\tau_{m} = (RC)_{mem}$ and we have previously set both $R = 1$ and $C = 1$, we need to set $\\tau_s \\geq 20$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = 20\n",
    "voff = 0\n",
    "a = 1\n",
    "\n",
    "i2 = neuron.add_current(a, voff, ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's now measure again the steady-state IV curve of the neuron. We can change the gain of the element and see how it changes the steady-state IV curve. Let's set the gain of the element so that the steady-state characteristic is back to a straight line:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "a = 2 # set the gain\n",
    "i2.update_a(a)\n",
    "\n",
    "V = np.arange(-3, 3, 0.1)\n",
    "I = neuron.IV_ss(V)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Neuron steady-state IV curve\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V, I)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The neuron's steady-state IV curve is again that of a positive resistor! But, importantly, the neuron does not act like a simple resistor anymore because of the *fast* nature of the negative resistor - in the fast timescale, the circuit is still a bistable switch. We can see this by plotting IV curves in *different* timescales. We define these IV curves as follows:\n",
    "\n",
    "$$\\mathcal{I}_x (V) = \\sum{I_j(V)}, \\quad for \\ all \\ \\tau_j \\leq \\tau_x$$\n",
    "\n",
    "This formula has a simple interpreation: if we want to plot the IV curve of the circuit in the timescale given by $\\tau_x$, then we just need to sum up the steady-state characteristics of all elements acting in this timescale or *faster*.\n",
    "\n",
    "To check the fast IV curve, we set the $\\tau_f = 0$ and call the `IV` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf = 0 # set the timescale\n",
    "V = np.arange(-3, 3, 0.1)\n",
    "I = neuron.IV(V, tf)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Fast IV curve\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V, I)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also observe what happens when we change `tf` in the previous code. From the equation, we expect that the IV curve plot will only change if we set `tf` to be larger than the timescale of the slow element. In this case, we would get the same output as if we directly measured the steady-state IV curve through `IV_ss` method."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The combination of the N-shaped fast IV curve, and the monotonically increasing slow IV curve means that the neuron will be spiking for a certain range of input currents:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trange = (0, 300) # time length of simulation\n",
    "I0 = -0.7\n",
    "i_app = lambda t: I0# Iapp(t)\n",
    "\n",
    "sol = neuron.simulate(trange, i_app) # call neuron's simulate functionality\n",
    "\n",
    "# Plot simulation\n",
    "plt.figure()\n",
    "plt.title(\"Spiking behavior\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Membrane voltage\")\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can play around with the value of $I_0$ to see how it affects the neuronal behavior. For a small enough input current, the neuron will stop spiking and be at rest. It is however still *excitable* and a short excitatory pulse can generate a spike. We can try different values of $I_0$ and pulse parameters to see when we get a spike at the output:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trange = (0, 300) # time length of simulation\n",
    "\n",
    "def pulse(t, t_start, t_length, mag):\n",
    "    I = (t>=t_start)*mag - (t>=(t_start+t_length))*mag\n",
    "    return I\n",
    "\n",
    "def i_app(t):\n",
    "    p1 = pulse(t, 150, 2, 0.3)\n",
    "    \n",
    "    pulses = [p1]\n",
    "    I = -1\n",
    "    for q in pulses:\n",
    "        I = I + q\n",
    "    return I\n",
    "\n",
    "sol = neuron.simulate(trange, i_app, max_step=1) # need to set max_step so that the pulse is not *missed*\n",
    "\n",
    "# Plot simulation\n",
    "plt.figure()\n",
    "plt.title(\"Spike output\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Membrane voltage\")\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Input current\")\n",
    "plt.plot(sol.t, i_app(sol.t))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creating rest-spike bistability\n",
    "\n",
    "We have so far seen how a spiking neuron can be created through the interconnection of a fast negative resistance element, and a slow positive resistance element. These elements replicate the fast excitatory nature of the sodium current, as well as the slower refractory mechanisms of the potassium activation and sodium inactivation.\n",
    "\n",
    "We will now see how introducing slower currents into the mix allows the neuron to modulate its spiking behavior so that it experiences more complex dynamical phenomena, such as bistability between resting and spiking, and bursting. These currents are usually dependent on calcium ions in biophyiscal neurons, and tend to act on timescales significantly slower than the processes controlling the excitability of the cell.\n",
    "\n",
    "We start with designing rest-spike bistability as it is a precursor to designing bursting behavior, very much how creating a bistable switch was the first step to generating spike outputs.\n",
    "\n",
    "Let us add another element to the circuit - this time it will be a slow *negative* resistance element! Importantly, this time we will need to change the $V_{off}$ parameter of the element to be non-zero, as the slower calcium-dependent currents activate in lower voltage ranges than the previously considered sodium and potassium currents."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "ts = 20\n",
    "voff = -0.9\n",
    "a = -1.5\n",
    "\n",
    "i3 = neuron.add_current(a, voff, ts)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's see how the addition of this element has changed the steady-state IV curve - remember, as the element is slow, it will not affect the fast IV curve, which we also plot for reference. You can change the element's gain and offset values to shape the steady-state IV curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Change the gain and the offset of the slow negative resistance element\n",
    "#i3.update_a(-1.5)\n",
    "#i3.update_voff(-0.9)\n",
    "\n",
    "tf = 0 # set the timescale\n",
    "V = np.arange(-3, 3, 0.1)\n",
    "If = neuron.IV(V, tf)\n",
    "\n",
    "Iss = neuron.IV_ss(V)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Fast IV curve\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V, If)\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Steady-state IV curve\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V, Iss)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Previously, when we had an \"N-shaped\" steady-state IV curve the neuron was bistable, and we were able to switch between the high \"up\" state and the low \"down\" state by injecting a series of current pulses.\n",
    "\n",
    "We now have the same picture, but with an important difference - the high \"up\" state now belongs to the negative resistance region of the fast IV curve, i.e. the high \"up\" state is now a spiking state!\n",
    "\n",
    "Let's test this behavior by again injecting a series of pulses and switching the system between the resting and spiking states. This time we will need to use longer pulses, as the circuit operates on the *slow* timescale, compared to the *fast* timescale of our simple bistable neuron."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trange = (0, 3000) # time length of simulation\n",
    "\n",
    "def pulse(t, t_start, t_length, mag):\n",
    "    I = (t>=t_start)*mag - (t>=(t_start+t_length))*mag\n",
    "    return I\n",
    "\n",
    "def i_app(t):\n",
    "    p1 = pulse(t, 800, 20, -1)\n",
    "    p2 = pulse(t, 1600, 20, 1)\n",
    "    p3 = pulse(t, 2400, 20, -1)\n",
    "    \n",
    "    # add more pulses\n",
    "    pulses = [p1, p2, p3]\n",
    "    I = -0.8 # this current needs to be in the *bistable* region\n",
    "    for q in pulses:\n",
    "        I = I + q\n",
    "    return I\n",
    "\n",
    "sol = neuron.simulate(trange, i_app, max_step=1) # need to set max_step so that the pulse is not *missed*\n",
    "\n",
    "# Plot simulation\n",
    "plt.figure()\n",
    "plt.title(\"Switching behavior\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Membrane voltage\")\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "\n",
    "plt.figure()\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Input current\")\n",
    "plt.plot(sol.t, i_app(sol.t))\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can play around with the gain and offset parameters of the slow negative resistance element and see how it affects the steady-state IV curve and correspondingly the switching behavior."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Bursting\n",
    "\n",
    "We will now turn the rest-spike bistable neuron into a burster. We can achieve this again by following the same procedure as when going from a bistable to a spiking neuron: we need to introduce a slower current that will \"automatically\" modulate the current so that when neuron is in the resting state, the current increases and leads to spiking, and when the neuron is in the spiking state, the current decreases and terminates the burst.\n",
    "\n",
    "We will call the timescale of this element *ultra-slow* as it needs to be significantly slower than the processes introduced so far. Let's define this element now:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tus = ts**2\n",
    "voff = -0.9\n",
    "a = 1.5\n",
    "\n",
    "i4 = neuron.add_current(a, voff, tus)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's measure the neuron's IV curves again. This time, there will be three IV curves of interest: fast, slow, and the steady-state (ultra-slow) IV curve."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "tf = 0 # set the timescale\n",
    "ts = 20\n",
    "\n",
    "V = np.arange(-3, 3, 0.1)\n",
    "\n",
    "If = neuron.IV(V, tf)\n",
    "Is = neuron.IV(V, ts)\n",
    "Iss = neuron.IV_ss(V)\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Fast IV curve\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V, If)\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Slow IV curve\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V, Is)\n",
    "plt.show()\n",
    "\n",
    "plt.figure()\n",
    "plt.title(\"Steady-state IV curve\")\n",
    "plt.xlabel(\"Voltage\")\n",
    "plt.ylabel(\"Current\")\n",
    "plt.plot(V, Iss)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When the input current is in the right range, the neuron will now burst!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "trange = (0, 2000) # time length of simulation\n",
    "I0 = -1\n",
    "i_app = lambda t: I0# Iapp(t)\n",
    "\n",
    "sol = neuron.simulate(trange, i_app) # call neuron's simulate functionality\n",
    "\n",
    "# Plot simulation\n",
    "plt.figure()\n",
    "plt.title(\"Bursting behavior\")\n",
    "plt.xlabel(\"Time\")\n",
    "plt.ylabel(\"Membrane voltage\")\n",
    "plt.plot(sol.t, sol.y[0])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can again try changing the applied current $I_0$ to observe the different spiking regimes of the neuron. You should notice that for low currents the neuron is at rest, but above a certain threshold the neuron starts bursting. Increasing the current further eventually terminates the bursts and leads to regular spiking, before again reaching a high resting state for a large enough input current."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Graphical interface with live plotting\n",
    "\n",
    "If you would like to explore further how changing the parameters modulates the behavior with a live neuron simulation, there is additional code that can be found at the following GitHub page [Circuit-Neuromodulation](https://github.com/lukaribar/Circuit-Neuromodulation).\n",
    "\n",
    "The script `gui.py` genereates a graphical interface for controlling the gain and offset values for the elements discussed here, and you can observe the different behaviors that can be achieved by exploring the different parameter combinations. Mainly, you should observe how changing each parameter affects the IV curves of the circuit, which in turn capture the main excitability properties of the circuit."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}