02-Building_a_neuron.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Building a neuron!\n",
    "\n",
    "In this section, we will implement one of the simplest mathematical model of a neuron: the **leaky integrate-and-fire (LIF) model**. The basic idea of LIF neuron was proposed in 1907 by Louis Édouard Lapicque, long before we understood the electrophysiology of a neuron (see a translation of [Lapicque’s paper](https://link.springer.com/article/10.1007/s00422-007-0190-0)). More details of the model can be found in the book [Theoretical neuroscience by Peter Dayan and Laurence F. Abbott](https://www.gatsby.ucl.ac.uk/~dayan/book/)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "\n",
       "        <iframe\n",
       "            width=\"560\"\n",
       "            height=\"315\"\n",
       "            src=\"https://www.youtube.com/embed/rSExvwCVRYg\"\n",
       "            frameborder=\"0\"\n",
       "            allowfullscreen\n",
       "            \n",
       "        ></iframe>\n",
       "        "
      ],
      "text/plain": [
       "<IPython.lib.display.IFrame at 0x7fcff127bc10>"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# @markdown\n",
    "from IPython.display import IFrame\n",
    "\n",
    "IFrame(width=\"560\", height=\"315\", src=\"https://www.youtube.com/embed/rSExvwCVRYg\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## LIF Neuron Model Description\n",
    "\n",
    "The subthreshold membrane potential dynamics of a Leaky Integrate-and-Fire (LIF) neuron is described by the following differential equation:\n",
    "\n",
    "\\begin{align}\n",
    "C_m \\frac{dV}{dt} = -g_L(V-E_L) + I \\quad \\text{(1)}\\\\\n",
    "\\end{align}\n",
    "\n",
    "where:\n",
    "- $C_m$ is the membrane capacitance,\n",
    "- $V$ is the membrane potential,\n",
    "- $g_L$ is the leak conductance ($g_L = \\frac{1}{R}$, with $R$ being the leak resistance),\n",
    "- $E_L$ is the resting potential, and\n",
    "- $I$ is the external input current.\n",
    "\n",
    "Dividing both sides of Equation (1) by $g_L$ gives:\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "\\tau_m \\frac{dV}{dt} = -(V - E_L) + \\frac{I}{g_L} \\quad \\text{(2)}\\\\\n",
    "\\end{align}\n",
    "\n",
    "where $\\tau_m$ is the membrane time constant and is defined as $\\tau_m = \\frac{C_m}{g_L}$. Note that dividing capacitance by conductance gives units of time!\n",
    "\n",
    "### Simulation of LIF Neuron Dynamics\n",
    "\n",
    "We will use Equation (2) to simulate the dynamics of the LIF neuron.\n",
    "\n",
    "If the external input $I$ is strong enough to make the membrane potential $V$ reach a certain threshold value $V_{\\text{th}}$, the neuron will generate a spike. In that case, the membrane potential $V$ is reset to a value $V_{\\text{reset}}$ (where $V_{\\text{reset}} < V_{\\text{th}}$), and the voltage is clamped to $V_{\\text{reset}}$ for a certain refractory period $\\tau_{\\text{ref}}$ ms, mimicking the refractoriness of the neuron during an action potential:\n",
    "\n",
    "\\begin{align}\n",
    "\\text{if}\\quad V(t_{\\text{sp}}) \\geq V_{\\text{th}}: \\quad V(t) = V_{\\text{reset}} \\quad \\text{for } t \\in (t_{\\text{sp}}, t_{\\text{sp}} + \\tau_{\\text{ref}})\\\\\n",
    "\\end{align}\n",
    "\n",
    "where $t_{\\text{sp}}$ is the spike time when the membrane potential $V(t)$ just exceeded $V_{\\text{th}}$.\n",
    "\n",
    "**Note:** In the video, the symbol $\\theta$ corresponds to the threshold voltage $V_{\\text{th}}$, and the symbol $\\Delta$ corresponds to the refractory time $\\tau_{\\text{ref}}$.\n",
    "\n",
    "The LIF model captures several important aspects of a neuron, including:\n",
    "- Spatial and temporal integration of synaptic inputs,\n",
    "- Spike generation upon reaching the threshold,\n",
    "- Refractoriness during the action potential,\n",
    "- Leaky membrane behavior.\n",
    "\n",
    "It is essential to keep in mind that the LIF model assumes that spatial and temporal integration of inputs is linear and that membrane potential dynamics close to the spike threshold are much slower in LIF neurons compared to real neurons.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Coding up a LIF neuron\n",
    "\n",
    "We now look at the Python code to calculate our equation for the LIF neuron and simulate the LIF neuron dynamics. We will use the \"Euler method\" to numerically integrate this equation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\tau_m\\frac{dV}{dt} = -(V-E_L) + \\frac{I}{g_L}\\,\n",
    "\\end{equation}\n",
    "\n",
    "where $V$ is the membrane potential, $g_L$ is the leak conductance, $E_L$ is the resting potential, $I$ is the external input current, and $\\tau_m$ is membrane time constant.\n",
    "\n",
    "The cell below initializes a dictionary that stores parameters of the LIF neuron model and the simulation scheme. You can use `pars=default_pars(T=simulation_time, dt=time_step)` to get the parameters. Note that, `simulation_time` and `time_step` have the unit `ms`. In addition, you can add the value to a new parameter by `pars['New_param'] = value`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# @title Plotting Functions\n",
    "\n",
    "def plot_volt_trace(pars, v, sp):\n",
    "  \"\"\"\n",
    "  Plot trajetory of membrane potential for a single neuron\n",
    "\n",
    "  Expects:\n",
    "  pars   : parameter dictionary\n",
    "  v      : volt trajetory\n",
    "  sp     : spike train\n",
    "\n",
    "  Returns:\n",
    "  figure of the membrane potential trajetory for a single neuron\n",
    "  \"\"\"\n",
    "\n",
    "  V_th = pars['V_th']\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  if sp.size:\n",
    "    sp_num = (sp / dt).astype(int) - 1\n",
    "    v[sp_num] += 20  # draw nicer spikes\n",
    "\n",
    "  plt.plot(pars['range_t'], v, 'b')\n",
    "  plt.axhline(V_th, 0, 1, color='k', ls='--')\n",
    "  plt.xlabel('Time (ms)')\n",
    "  plt.ylabel('V (mV)')\n",
    "  plt.legend(['Membrane\\npotential', r'Threshold V$_{\\mathrm{th}}$'],\n",
    "             loc=[1.05, 0.75])\n",
    "  plt.ylim([-80, -40])\n",
    "  plt.show()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'V_th': -55.0, 'V_reset': -75.0, 'tau_m': 10.0, 'g_L': 10.0, 'V_init': -75.0, 'E_L': -75.0, 'tref': 2.0, 'T': 400.0, 'dt': 0.1, 'range_t': array([0.000e+00, 1.000e-01, 2.000e-01, ..., 3.997e+02, 3.998e+02,\n",
      "       3.999e+02])}\n"
     ]
    }
   ],
   "source": [
    "def default_pars(**kwargs):\n",
    "  pars = {}\n",
    "\n",
    "  # typical neuron parameters#\n",
    "  pars['V_th'] = -55.     # spike threshold [mV]\n",
    "  pars['V_reset'] = -75.  # reset potential [mV]\n",
    "  pars['tau_m'] = 10.     # membrane time constant [ms]\n",
    "  pars['g_L'] = 10.       # leak conductance [nS]\n",
    "  pars['V_init'] = -75.   # initial potential [mV]\n",
    "  pars['E_L'] = -75.      # leak reversal potential [mV]\n",
    "  pars['tref'] = 2.       # refractory time (ms)\n",
    "\n",
    "  # simulation parameters #\n",
    "  pars['T'] = 400.  # Total duration of simulation [ms]\n",
    "  pars['dt'] = .1   # Simulation time step [ms]\n",
    "\n",
    "  # external parameters if any #\n",
    "  for k in kwargs:\n",
    "    pars[k] = kwargs[k]\n",
    "\n",
    "  pars['range_t'] = np.arange(0, pars['T'], pars['dt'])  # Vector of discretized time points [ms]\n",
    "\n",
    "  return pars\n",
    "\n",
    "\n",
    "pars = default_pars()\n",
    "print(pars)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjAAAAEiCAYAAADnB+hiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAACIZklEQVR4nO2dd3xT5ffHP0kz2iSd6aBl7z3LUNmgsocLRURRUFBx4Ub9On+Kigv1q1++KqI4vjgAQQSRoaAgsmTILqNA6Uyaps3O8/vj8twmbVo60twbct6vV143vbm597m3N+f53POc5xwFY4yBIAiCIAgijFBK3QCCIAiCIIjaQgKGIAiCIIiwgwQMQRAEQRBhBwkYgiAIgiDCDpXUDSAIgiCIhoIxBq/XC4/HA4/HA6/XK/7t+973b8YYGGPiOv43f1VEoVAEfCmVSkRFRSEqKkp8r1QqxZdKpRLf+24XLvBr5HK54PF44Ha74Xa7xevo++L/B76seF0rolAoEB0dDYPBAKPRCIVCUXkbmoVEEARBhBLeiQHlnX+gbfjS9xWoM2SMwe12i0KEd6R8nUKh8BMIFcVExc94m/g6LkZ4e/myqjb6ttVXHPmKJK/X69fZ8/e+x61K6AR6X7HtVV1T3uaKAqKi4PB4PAHX+15nj8eDqKgoqFQq8cX/5tfM9/oFel9dW202GwoLC6FSqZCenl5pOxIwBEEQRMgoKyvD6dOnxc6ouidwvvT1aFR8H0iEqFQqqNVqsUMNF6+Gr+ipyjtU1fuaejX4da0oICoKuKioKL/r63ud+XXl2zQkXq8XJ0+eRGJiIhITE/3PhwQMQRAEEQrcbjeysrLQqFEjxMXFSd0cIkwoLS1Ffn4+WrRo4bc+PGQpQRAEEfa4XC6o1WoSL0StUKvVcLvdldaTgCEIgiBCAmMsbIZzCPlDdxJBEAQREniQKkHUBsZYwPuGBAxBEAQREqrqiMKFzZs3w2Kx+K07cOAAfv/99wY75mOPPYZ//vmnwfYfDni93oCeOxIwBEEQREioqiMKBw4fPoxBgwZh/vz5fuunTZuG6667LmCMRjB44403kJeX1yD7Dhf4dO2KyPJO+ueff3DixIlK6xlj2LVrF7Zu3QqXyyVBywiCIIi6UlVHFA6sXLkSALBx40ZxncViwc6dO5Gbm4ujR48G/Zh8WrVarQ76vsOJsPHArFq1Cl27dsWdd97pt/7AgQMYMGAAMjMzccUVV6BLly7YvHmzRK0kCIIgaks4e2D+/PNPAMCePXvEJHxbtmwR860cOXKk0nccDgfOnTtX6YH7gw8+QEFBAbxeLw4ePCh6b4qKimA2m8Xt+L5VKhXKysqQlZVVydPDE84BQHFxMY4dO1apHVarFYWFhX7r1qxZg5ycHDDG8NdffwX8ntPpxLlz5+BwOKq+MCEgLDww+fn5mD59OnQ6HTwej7i+pKQEI0eOhNfrxe7du5GVlYVevXph/PjxlcYjCYIgCHni9XrD1gOzZ88eNG/eHFarVezsV69ejUaNGqFRo0Y4dOiQ3/bvvfce0tPT0bhxYyQmJuLVV18FIIiJe++9Fxs2bMCoUaPQqVMnjBs3DtOmTYPRaETjxo2xYcMGv3299tprSElJQevWrdG+fXusWbNG/Oyaa67B22+/jTvvvBNpaWlo27YtvvnmGwBATk4OJk+ejMTERCQnJ6Nfv37IysoCADzwwAN49913ce2116Jv377o2LEj1q5dK+538eLFaNasGRo3boyEhATMnTs3+Be1hlTpuWMywev1smuuuYZ17dqV3XnnnWzIkCHiZ++++y7T6/WssLBQXGe1WpnBYGDvvPOOFM0lCIIgasm5c+dYUVGR1M2oNfn5+UyhULBXXnmFAWBffvkl83q9rGXLluzee+9lV111FbvtttvE7X/88UemVqvZf//7X3b48GE2b948plAomNVqZSaTiQFgHTp0YO3bt2cPPfQQA8BatWrFlixZwtq3b88mTpzIGGPM4/EwAKxly5bs66+/Zrt372Z33nkni42NZXl5eYwxxvr06cMAsCuvvJL98ssvrGfPnuypp55iZWVlrEOHDqxNmzZs9erVbOvWrQwA+/bbbxljjHXp0oWp1WrWvn179scff7DLLruM3X777Ywxxv744w+mVCrZO++8ww4fPsz+/e9/MwDszJkzob3wFzh//jwrKCiotF42xRyXLFmC5cuXY9u2bVi8eLHfZ2vXrsU111yDpKQkcZ1er8cVV1yBv/76q9K+Dh48iEOHDkGpVEKr1SIhIQE6nQ5du3YN6wh4QHApFhcXo7CwEMXFxSgtLUVxcTFMJhMKCwtRUlICh8MBp9MJp9MJl8uFsrIylJaWwmazwel0irUsfOGpo1UqFTQajZiGW61WQ61WQ6fTISkpCXFxcYiNjUV8fDz0ej0SEhIQHx+P6OhoREdHQ6/XIz4+/pIds3W73TCbzbBarSgtLYXFYhGvrc1mg91uh9VqRUlJCcrKysSX0+mEw+GA3W6Hy+Xyq9XimwIcKE+hzq+777XVarVQq9UwGAyIj49HfHw84uLiEBcXJ75PTU1FfHx82N7rJSUlKCoqQmlpqfgqKytDSUkJSkpKxOvL3/Nrarfb4XA44HK54HQ6/e5xhUIh3tsajQYxMTGIjY0VX77XLyEhAQkJCeL7xMTES+J+5sMZJpMJRUVFyM3NFe9fu90u3qsOh0O8p/m9ylPWd+vWDa+//nqd2+B2u6FSyabbqTHr168HYwzXXHMN/vvf/+L3339Hu3btcOLECdxwww1Yvnw5Nm3aJG7/xhtv4L777sPNN9+M//znP3j99ddx9dVXQ6/Xi16o0tJSbNu2Dfv27cNbb72F7du3w2g04tixY6IHhf+Gv/jiC1x++eUAhOGn9evX4+uvv8Z9990Hj8eDfv36YfXq1VCr1Rg9ejQGDhyIpUuX4sSJE9i/fz/atGkjxui0bNkSAMRhoRUrVqB9+/bo1asXjh8/DgB48803ceutt2LmzJlYtGgRXnjhBfTr1w+NGjUKyfWuiMfjgVarrbReFndSdnY27rvvPsyZMwd9+/atJGAOHjyIadOmVfpeamoqzp07V2n90qVL8dxzz1Va73a7MWfOHOzfvx8xMTFISEhAUlKS2CHHxMTAYDAgMTFRNGZJSUlISkqCXq8P2g/P6/XCZrOhpKQEFosFZWVlsFgssFgssFqtyM3NRW5uLs6fP4/CwkLxM5PJhJycHNjt9mr3r1AoREPNjbVer0dMTAy0Wq1fDQtekMzj8cDhcMDtdovCx+12w+VyiSLIbDaLY78XIzo6GgkJCTAajTAYDNDr9UhKSkJycrLYMaSmpsJoNEKv14sdCO84YmJigt4BO51O5Ofno6ioSOz8CgsLUVhYKHaEVqsVJpMJFosFxcXFKCkpETtRq9UqjlvXlJiYGMTExECj0UCr1SI6OloUh/zlW2cEKA/cs9lsyM3NFYVRWVmZ2Fk7nc5qj6vRaJCamoqUlBSkpqYiPT0daWlpSEtLg06nQ0JCApKTk0XXckJCAgwGQ9DiExhjcDgconjmIoSL75ycHJw/f15cnj9/HkVFReL/oiZotVoYDAbExMRApVIhOjpaFHgajcavTovH44HdbheFvd1uF39/NpvtosfS6XQwGAyIjY0Vr6nRaERSUhJ0Oh1SUlKQnJws3uvx8fFITEwUxVAwritjDE6nE2VlZbBarbBYLMjPz4fJZBL/5ufEH2pycnKQn5+PvLw85OfnV7v/qKgo6HQ6aLVa0V743qtRUVEoKyur1zmEaxDvzz//jLS0NLRr1w6DBg3CL7/8Aq1Wi7S0NAwYMAAnTpzAggULUFJSgtjYWOzduxcpKSlo0aIFYmJi8NJLL2H69OkAINYTeuaZZ5CRkYGTJ08CgHhd2rRpg6ysLHHKOX+g5ERFRaFLly6i2ACACRMmiNu89NJLAIBbbrkF/fr1Q5s2bQAINgEQbDMAmM1mXHvttWjfvj0AiBWlAWDv3r3o0qULWrVqBa/Xi8ceewx33323ZP+7qmKnJBcwTqcTU6ZMgdFoxNNPPx2wc9BoNJU8BoCgIAM9GQXalv8AAcBut8NsNuPAgQMwm80oKSkJ+J2KqNVqaLVaaDQa6HQ68elYq9X6FQ3jRbZ4R+NyucQOkHdCFyMqKgqpqalITU1FbGws0tPT0bFjRzRq1Ajp6elITk4WvSDx8fFISkpCYmIi4uLixEqgwcbr9YpPwmazGaWlpTCbzSguLobdbofdbhc9QvwpuqioSPRW7Nu3D0VFRbBYLBcNCouKioJerxcFGO+kuEcoUPl530q0vJPinZbVaq1Rx8g7d+7diI2NRVpaGvR6PWJjY8X/iV6vF9dxY89fvKOLjo5usIBFl8sFi8UCs9ksdlzFxcUoLi5Gbm4u8vLykJeXh4KCAuTk5GD//v3Iy8urdvaeQqEQxSMXAWq1WrzHfYu7KRQK0XPkdDphs9nEjpU/vbOLlFlTKpVITU1FRkYG0tPT0bVrVyQlJSEjIwNGoxE6nU68zjqdTvT+GQwGGAyGoHlFPB6Pn2A1m83idTWbzTCZTKKdKCkpQV5eHk6dOoUdO3bAbDZftFPn11Wv14vXldsR34KDvC38HnY4HHA4HLDZbKLX72LXFBBsHX8YSEtLQ/v27dG/f380btwYjRs3FoVrWloa4uPjRTumVqtDUpgvHIN4N23ahEGDBkGhUGDYsGH49NNPkZeXhxtvvBFRUVHo06cPvF4vtm7diquuugoWiwV79+7F/PnzMXnyZPFe5Q8mTqcTaWlpAIQHcUCIV0lISEDTpk1hs9mQk5ODjIwMaDQav37R4/Hg0KFDGDRoULVtjo6OFkULAKSkpAAo97yYzWYMGTJE/DwqKkoU8xaLBXv27MGzzz6LW2+9VRQ9Uv3/ZCFgxo8fjz179sBms8HhcODdd9+F3W4XZxP5VppUKBRo2rQp1q1bh9TU1IDz4M+ePYuePXtWWt+xY0dMmDBB9CqYTCax833nnXcqbc8YQ1lZGWw2m/gEXlxcDIvFgoKCAphMJvEJkg/PcHcrd1tzNytjDCqVyu9JkLv8+dMiN8jcfc09EHFxcTAYDEhJSYHRaJTdEIBSqRTPIz09vV77KisrQ15ennhteefr2yFbrVaxc+SeB/7yrdLKjTp/6uZDBXzoRaPRwGAwICkpSXxS5h1hYmIiUlJSoNfrAwoOpxPwsQEh4+efgfHjgZgYYOJEYNEi/8/VajWMRiOMRmON9+n1esUhAz6MwD1QvtefDx1w8c3vcX6t+Yt7j7j3kos3fn/ze53/ze9zo9EoCkE5dGZRUVEBK93WFK/Xi4KCAtF75DusazabxZklpaWl4v3LhTX3dPIOyvce1mq10Gq1oqg2GAyIjo4WbQe/lklJSTAYDKLAawjvJScjIwMTJkzABx98UKfvh6OAKSwsRFZWljgzdvjw4QAAk8mEa665BgDQoUMHGAwG/Pzzz7j66qsxZswYFBcXY9iwYVCr1SgqKsIrr7yCAwcO4JNPPgEAxMbGAoAoZLKzs9GxY0dkZGQAAE6fPo2MjAwolUqsW7cO7du3h9vtxgsvvIDs7GzcfPPNYhsDCdtu3brhq6++QnZ2Npo2bYqsrCxERUVh586d6NGjB1wul19NqtTUVNEbNGbMGOzbtw/Dhw9HdHQ0LBYL3nnnHaxevRpbt24N5uWtEVUOPTZ8+E05x48fZzt37mR79uxhR48eZS6Xi5WVlbEtW7awdevWsR9//JF9//337Oqrr2bt27dnX375JbPb7eyhhx5inTt39tuXxWJhGo2GLV68uNbtaNu2LZsxY0awTou4xJgzh7GUFMa0WsbS0qRpw+rVjAHCa/BgadrAGGP5+YxZLIw5HIx5PNK1g5AHLVu2ZFOmTKnz9w8dOsTcbncQW9Tw7N+/nwFga9asEdf16tWLtWzZkrlcLnHd1VdfzTIzMxljQl/Xq1cvBoAlJiYyhULBWrRowbZs2cJsNhszGAxs6dKljDFhAotarWafffYZY4wxm83GALD//e9/jDHGYmNjmUqlYgAYAJaamspWrFghHjczM5O9/PLLldptsVhY+/btmU6nY+3atWNarZY1b96cDRw4kDHGWFJSEvvPf/4jbr9o0SKmUCiYx+Nh586dY/3792cAWEJCAlMqlSwjI8PvGoSSw4cP+11rTkg9MK1ataq0TqVSoX///n7r1q9fD6fTicmTJwMAJk2ahLfeegvfffedmPHw4YcfBgBceeWVtW6HRqNBUVFRHc6AiAQKCgAeLlALB0dQ8X2gMpmkaQMAdO4McOfn9dcDF2ILiQhFp9PVKGYoEOxCbFe4eWDatm2Ljz/+WPS8AMCyZcsAwM8r8H//93/Ys2cPAKGv27FjB/7++2/Rk9KrVy/x3H/66Sf06NEDAMRhqaZNmwIQhn7efPNNZGZmAhCmb6ekpGDXrl1wu90YMGCAX0DrSy+9JMax+BIbG4v9+/dj48aNyMvLw7Bhw6DVasVg45UrV6JTp07i9uPGjcOHH34IpVKJ9PR0bN68GQcOHMCJEyeQkpKCPn36SBIDw++bQMdWMFaDQdUQc8899+DgwYN+GQ8feughvP322xgwYAAKCwtx6NAhvPHGG3jooYdqvf8+ffogJSUFq1evDmaziUuEGTOAjz8W3vfpA2zfHvo2rF4NjBkjvG/ZEriQuiHkpKcD588L78eNA374QZp2EPKgPrbT6/XiyJEj6NChQwO0jLhUcbvdyMrKQrt27Sp9JnkQbyBmzZoljsVx3nrrLYwbNw7ff/89tFotvvnmG3Tu3LlO+9doNJJnFiTki+9Q64Vh6pDjG8JwkckjDYrvQ08dH7yJS4j62M5wnUJNSEt1940s76Zu3bqhW7duldYPGzYMw4YNq/f++UwhggiE7+QWqQSMr5fdahWGlKSI6fa1GzTqStTHdobrFGpCWqq7b8JrMDJIKJXKGk1HJCIT33xJ9ZxsVWcqhgmUlEjTDl8BQ05Loj62s6r4F5fLP+aLIHwJOw9MQ8MYC7tAMiJ06PXl7+PjpWlDxdszNxfwmfEYMny9URfJnUdEAPWxnYE6IsaACylGEBMDJCQASUmC5zM+XlhnMAjrExKE34DRKGyj1wufa7VCqgOdTnhFRwvrVCphCFSpBLxewOMRxJLLJfztk9idkDHVeWAiUsB4vV4aiyWqxHfYSKeTpg0Vb896JkCtM74Cpgb5F4lLnPrYTq/XWyk/DS+s7PUCpaXC6+zZ+rayZlzKXp+azvb6/fff8fnnn+PDDz8MQavqRnXnEpG9OLuQopkgAuHrdbmQJDPkVOwjpBIPBkP5e6mGsQj5UB/bGch7o1aXe0bsdiHOymwGLBbhZbMJMWBFRcLfxcVAYaGQWqC0VBD2DofgHSwrE152u7DO7Rb2zYmKEo6n0Qjvvd7Kns6G4tSpU2CMoUWLFkHd7+bNm3HFFVf4eSieffZZvPjii8jJyRGT5FXFjh078N1338lawHg8HhpC8oWCyYjqSEgofy+Vm7lilvwalgcKOr4eKKnaQMiH+thOt9td5XfVauEVGws0b16fFsqT5557Drm5uUFN3WG32zFo0CCsWLEC48ePF9ffdtttcDqdSE5Ovug+wqEvDOS540RkIIjD4QhY2ZIgAP8hJKmCeCsKmMJCadoRE1P+3u0WnlqJyKU+tpMxJvvOsip2794NNx/vqgXsQtkNjUYDk8lUqwDokydPinXzeP0k3/3yOkcV99uqVSu88sorNbrWLpdL9pXWq0piB0SogLHb7WJxKoKoiO8QUh3L49Sbin2EVN4PXwEDCK5+InKpj+30eDxhOXkiPz8fmZmZePfdd3HdddehZ8+emDdvnp+gOXz4MG6++WZ06NABs2bNgtlsBgDMmTMHOp0On332GXbs2AGtVguj0SgKkx07dmDEiBFIT09Hnz598MUXX4j7HDduHBYsWIBXXnkFCQkJaNasGXbv3g2n04mkpCSxDtrdd98NpVIpZqdfunQp/vWvf4n7cblceOedd9C9e3e0aNECM2fOFDPRh4uAqeq+Cb+7KQiEwz+NkA5fz6vvjKRQUrGPkCqIt2IeHPLARDb1sZ3hGnvodDrBGMOcOXNgtVrRu3dvPPXUU/jPf/4DADhx4gT69u2LXbt24eqrr8ZXX32FBx54AIAgLubPn4++ffsiOTkZ7733Hj788ENotVpkZWVh8ODBSElJwRtvvIFBgwbh9ttvR05ODgChyOj8+fMxb948zJ8/HwkJCfjPf/4DjUaDxYsX4+mnn4ZGo8HYsWPx8ccfY9q0aQCA3377Tcxi7/V6cf311+PJJ5/Etddei3/9619YuHAhVq5cKX4ud68YzUKqgNPp9CszThC++HbaUk2jrjhlurhYmnZUDGK+lGduEBenPrYzHOItAsGrlL/wwgt45plnAAjxPP/73/9w77334tVXX0ViYiJ27doFnU6HZs2a4ZlnnsHChQvRrl07tGvXDufPn0d+fj7uuusucb/z58/H5Zdfjs8//xw7d+7EwoULkZycjPgLRsftdiM/Px/r1q3DlVdeiSNHjuDIkSMAIMa8zJs3DwMGDMAdd9wh7tdisYiVrjdu3IgffvhB3MeZM2cAAK1btxaPIfe+sDrhSx4YgqiAb+CqFLlXAh1XqiRyUhWzJORJfWxnOBZyBICYmBgolUq/YsRdu3YVy9388ccfmDRpEnQXDMfYsWNht9uRm5vrtw97hamEW7duRadOnTBmzBj069cP6enp2Lp1q7gfs9mMoUOHigWLvV6vXxxMVfstKSkRBcx3332Hrl27ivvg119/wbUcDuEU1QmYiPTA2Gw2xFQc3CeIC8ghzkOtLk/ABUjXJt8ZWYA05QwI+VAf2xmuAkahUCA6OtqvCrfdbhc9F74eDwDidok+AXSBhEZubi4+/PBD3HbbbTh06BDatm3r97nJZMLYsWPFv71eL1wVDEGg/fp6ulwuF1JSUsTPki5Mq/RcmF9eVlYW1gIm/O6meuL1emGxWJBQ0TITxAV8851IlX9FofAfvpLKA1MxBoYETORSX9sZrjEwgBCPkp+fD8YYzp8/j08//RRXXHEFAGDgwIFYvnw57HY7vF4v3n33XWRmZvqJGq1W6yeAAOCyyy7D6NGj8d5774niZe3atZg6dSoAQST5io/k5GRxCKi6/arVatFT06ZNG+zZswelpaUAhCHA6OhoHD58GIAgYOQ+hERBvD5YrVYwxsRxRoKoyIUJBACEhFlS4dtPSCWkKv5MKIF15FJf2xmuHhhA8MLMnTsXBoMBGRkZsFgseP755wEAjz76KI4cOYLWrVujTZs2+Oyzz/Daa6/5fd/lclU69+effx5bt25Fs2bNcNVVV6FTp04YPXo0evbsCQCIi4sThQcANG/eHNnZ2X7FNAPtV6/Xi0Jx2rRp8Hq9yMzMxJQpU9CxY0e43W5xtpPL5RKHrOQKDSH5wKe3kYAhqsJXtBQUAE2bStMOX0+9VHWIKibyC8MYTCJI1Nd2hrOAYYxh9uzZaNeuHeLj43HdddeJcSTdunXDnj17sGjRIpSWlmLmzJno0KGD3/f79+8vCh5O165dcezYMXz//fc4ceIEUlJSMHr0aLRs2RIA8Oqrr2Lo0KHi9sOHD8f06dP9OvPbb7/db5gJAF5//XVxiCgtLQ379+/HkiVLkJubi48//hgtW7bEkiVLAAD33ntvpWEpuUGlBHwoKCgAAHEOPUFUxGIpf+/rjQk1vv1EHXJoBQVfARMTQ0NIkUwwbGe4DiExxtCjRw9Mnz494OetW7fGSy+9VOX3e/fujd69e1daHxsbi9tuuy3gd3xnLAFAy5YtsXDhQr91jz32WKXvVSwf0LhxYzz++ON+61588UUAwjCW3KEYGB9MFx6vScAQVeErYKRMn+97i0r1kOTbBop7j2zqaztrk4VWblSXzp6QjogTMPwpIolqqRNV4Ot1kVLA+GYB9i1KF0p8p3NT5oHIpr62M5yDeLVaLc1clSERN4TEx3ETpcoRT8geX9GSlyddO3z7CakeXinmheBEsu3cu3cvGjVqJHUzIhYK4r1A2YWc7HqpcsQTssd3CEmqIoqAfwyMHLzvUnmBCHkQybazSZMmUjchoqnKexdxQ0i5ublQq9WIkyrFKiF7fPPASBnE6/ugKwcBI/PJCkQDEwzbGa5DSIQ8iUgBk5qaGrbT+YiGxyf1gmQ1iAB/ASMHu+97XYjIIxi2M5wDeQn5EXG9eE5ODo1lEtXim3NFyhgYuQkYqaZyE/KAbCchBQqFokrhG3ECJi8vD+np6VI3g5Axvmn7pRQwPlnEIQeHoRzaQEgH2U5CbkScScrPz0dycrLUzSBkjG+sh288TKjxTbchh9lAMq/5RjQw9bWd1T1JE0RVkAfmAowx5OXlITU1VeqmEDLGd7aNlHEfvrWQ5OD9oDpIkUswbCcF8BJ1gQTMBYqLi+F0OknAEDXGd0p1qPGdRi0H8UCzkCKXYNlO8sAQtYUEzAXyLgQ0VKwVQRBVYbUCPsVfQ4pWW/7eZpOmDb7IoQ2ENATDdobzENKRI0dw/vx5eDweMMbgdDrhcrn8zmf37t246aabQnqOLpcL48aNw/Hjx+v0/UOHDmH48OFwVxOh/+233+KTTz4J+JnNZsOJEycCfnbw4EGcO3euTu3yhQTMBSwXHqepEjVRHRWHa+x2adrhC4kHQkqCYTuVSiW8Uj0N1IPz58+jffv2SE9Ph0qlglKphFarhUajQXR0NFJTU9G7d29kZ2fjf//7n3itQsVPP/2ErVu31um7x44dw4YNG8QkhYFYvnw5Nm3aFPCzTz75BH379q30fzWZTBg0aBAWLFhQp3b5Up2AkYFjOnQUX0jqQQKGqI6KNX8sFkCnk6YtHAofIKQkGLYzXD0wjRo1wt9//439+/fDarVi0aJFKCwsxBNPPIHS0lKYTCbExcUh4ULQmtM3D0MDo1arodPp6nxMz4WAP1U1Y9SFhYVo2bJlwM/atWuHgoICHDp0CJ06dRLXf/TRR7BYLLjvvvvq1C5fSMBcgCvj2NhYiVtCyJlAAkbq9BdySOMvh0BiQhqCYTvDVcAAQLdu3dCtWzcAwN9//41t27bhjjvu8Ntm7dq1AFDjoo8ulwsqlapWwc25ublYtmwZ9u7di7fffhsajQZ2ux26Oj5hceGjrqZSq9VqrXL2Wf/+/aHVarFu3To/AbNkyRJcc801aNy4cZ3a5Ut1nruIMkn8R0hlBIjqqPhblsPwjRwCaOUQSExIQzBsZ1RUlPjEH84YjUaYTKZK63kna7FYMHXqVLRp0wa33HILbBcMyKZNmzB58mRs2rQJw4cPh0ajwdSpU8XvL1++HF27doVOp0OXLl3wzTffiJ/t3bsXQ4YMQXp6OmbPno0PPvgAR44cEY+rVqvx8ssvo2PHjhgwYAAOHjwofrekpATPPPMMmjdvji5duuCnn34SP+OxL9wDwxjDf/7zH3Ts2BGxsbEYM2YMzp07B6NvTgcfdDodBg8e7LfPAwcOYO/evZgyZUrtLmwVUAzMBbgbNMF3fipBVKDiA5QcUujLIQsuOS4jl2DYzuqepLVaLQwGA5o0aYJu3bqhf//+GD16NK677jrRswEIgaGffvopfvjhB6xfvx5//PEHdu7ciX379uH48eN+HZ3b7RaDbj0eD5xOJ0pLS8Wq2nXFaDSiqKio0np+7OHDh2P79u0YO3YsvvjiC6xatQoAcPToUXzzzTcYO3YsevfujdmzZ2P9+vUAgF9//RXXX389rrnmGnz77bcYNmwYpk2bJgYLX3vttTh79iy++uor5Ofno0ePHmJJB8YYnn76abz88su48sorcfr0acybN09s15QpU/D2229j0qRJSE9Px8SJE3H27FkAgvhRKBSiF+idd97Bo48+iltuuQUfffQRGGPIyspCUlJSlddj1KhR+PXXX2G/ECy4ePFiJCcnY+TIkfW6zpzq7puIeqbiP0LywBDVIUcPjFQPrr52I0y9/0QQCIbtrOpJ2u12w+l0igKDd66csWPHiu83bNiA2bNnV7l/Xw9Pjx49cODAgYDb1mcoy2g0wmKxwOv1+tWF4vuMiorC9u3bERsbi40bNyIrKwuAEG/i8XiwdOlSjB49GuvXrxdnD7388su488478dxzz2HFihVYs2YNunTpArVajYKCAhw/fhwzZsxAixYtEBcXh927d/sd9/jx49i6dSsyMzMRHx+PX3/9FQCwc+dOrFy5EqtWrcKYMWNQVFSEVq1aYdmyZZg9ezbcbje0F6Y7Msbw9ttv491338Vtt90GAJg4cSKSkpJgMBiqvB4jR47EQw89hC1btmDo0KFYsmQJbr75ZnFY6vjx4zhz5gwGDx4MAHjxxRfRqFEj3HnnnTW63lFRUSRgAGEsT6PRVDveRxC+05cBeQzfSBXE65uJWA6eKEIagmE7qxpCUqlUcDqdsNvtKCoqgtlshsVigcVigc1mQ+/evcVt27dvj6lTp8JkMqG0tBRlZWVwOBxwOp2V4kkqxpZERUVBrVZDo9FUEh+1wWAwgDGGkpISv6BmLmDee+89MVZIqVSixOdH1KJFC4wePRqA4KkZPnw4AOCvv/5Cly5d0KNHD2RlZWH27Nl48sknoVAokJycjLlz5+L999/HRx99hDZt2mDJkiXo16+feNz7778fmZmZ4nnzY/7555+Ii4vDqFGjAABJSUnIzMxEdnY2AGEaNBcw2dnZOHXqFIYNGya2V6vVIioqqto4nfbt26N58+ZYvXo1PB4PcnJyRAEEAPv27cP69etFAeN7vRhjF/0/KBQKEjCAEDRF4oW4GBXj4XxrI0lFCCc2+OFbC8rhELwwNCMq8giG7awuBkatVkOtViM2NhbNmzevch9XXnklrrzyyhodb9++fXVq58XgQbp2u91PwPBO1ldwqdVqv+GmQJ212+2GyWTC4sWLcf/992P27Nl+QzYKhQL/93//hxdeeAG7d+/GrFmzcPvtt/t5l3yPqdFoUFhYCAAoLS1FdHS033GtVitSLhRaKysrQ/SFGiFcpFScUq1QKOCoxggqFAqMGjUKK1aswPnz59G5c2f07NlT/Pytt95CdnY2zpw5g3vvvRcAsHHjRnz11Vc4d+4c5syZg7vuuqva/VMMDACHwyH+swiiKirGeshBwEhVC6liHio5xOIQoScYtvNSCeLl3pWKsTR8Ro/Wx4WbkZHhNyQWqCNWqVS44oorMG7cODz++ONISkqCx+PBggULMGvWLGzbtg1//PEHoqKi0Lt3b8yaNQsHDx6EzWYLeMz09HScP38ebrcbw4YNQ15eHn788UcAwLp167Bjxw5RBJaVlUGj0QAAmjRpgq5du+Kxxx4TRcyyZctQWloqDiFWxQ033ICsrCwsXboUt912m5/H5qGHHsKYMWOwbNky8bgejwfr16/Htm3b8Morr1S7b4qBuUBpaWmdp5sRkUNFASMHmytV3HlOjv/fbnflGCHi0icYtjOcp1H7woc9KsYDFRYWIi0tzc9T1aFDB3z//fd+3w3E/PnzMWHCBKSnp6N169Y4efIk7Ha7GAT8f//3f2jevDk0Gg2OHTuG6667DjqdThRHTZs2FffVsWNHuN1uHD9+HL169cKkSZMwfvx4NG3aFKdPn8b06dPRo0cPAIKQ4PEtCoUC77//PsaMGYNGjRohKSkJp06dQnx8PHIqGoIKDB48GGlpaSguLsa0adMueg2HDBkChUKBhISEi4paygNzAbvdTh4Y4qJUjFeTg4CpJoauQalYjdvprDxLi7j0CYbtDNdMvBW54oorcOzYMaSnp/utnzx5Mi6//HK/dXfccYeYC2XIkCFVdtaXX345srKy8NNPP+HMmTNo0qQJRowYgbi4OIwbNw79+/fHli1b4PV60atXL0yYMAGA4OHZsmWLGP8CAH379sWzzz6LRo0aQaFQ4KuvvsJ1112H3bt3Y8CAAWIMDgDceuutGDhwoPj3wIEDcerUKSxbtgwWi0Wc7l1dEC8geNfmz5+PkpIScXiKE0iA+CbOu1genKioKLiqCESMOAFT0yRDROSi1/v/LZXN9T3uBS9vyKmYFV0OAc1E6AmG7VSpVNXW3AkXFApFwMy0cXFxYrI7Trt27dCuXTsAgjemQ4cOVe7XYDDghhtuCHi8UaNGiYG4FT/r37+/3zq1Wo3nnntO/FupVGLSpEmYNGlSpe+3bNmy0rkkJiZWStJXE2655ZaA61u2bImtW7di1apVouenNlAQ7wXKyspIwBAXpWK2dKm83lZr+XupgngrTiGXgzeKCD3BsJ2XigeGqB3dunXDq6++ir1796JZs2YYPnw49D5PiXPnzq32+xQDcwGahUTUBLmESfnmypKqoGTFISQSMJFJMGznpRIDQ9Se6maPzZw5s9rvKpXKKofeImoWEhB4GhtB+FJRwEg1bdh3koNUM6FCXFiXkDH1tZ3kgSHqQnVDj7LwwPzrX//Ctm3bxPnmjDG43W5MmDABjz32GABhrvzChQvx3nvvweFwYMqUKXj00UdrVVyM1D9REyrGq0k1hflCKgcA0sWe1DPrOnGJEAzbealMoyZCi+wz8RYUFGDHjh2YNGkSNBoNGGPQarVi0BJjDDfeeCN+/PFHzJ49G/Hx8fj3v/+NjRs34tdffyWvChFUKnpgpAqgzc0tfy9V7GNFDwwlsSPqCnlgiLog+xgYj8eD7t2748MPPwz4+Zo1a/D999/jt99+E6d8TZo0CR06dMDPP/9c46JRFWtlEEQgKjr1KpYWCBW+uaOkch5WzF9FzwqRSTBsJ58uW580/kTkIfs8MLm5udBqtXjwwQexbds2uFwuzJgxAzNnzoRSqcSyZcswePBgv/nq7du3R9++fbF8+fJKAubgwYM4dOgQlEoltFotEhIS0K9fP3oCIGqEwSB4YXhGbalSB1UMoJWCivWPKAY+MgmG7VQoFNBqtbBarVRQl6gxDofDL2+ML7IQMDk5Odi+fTuOHTuGG2+8EcXFxWKlzPvuuw9//PEHxo8fX+l7zZs3x6lTpyqtX7p0qd88eEBI80wChqgJCgXQqBFwoYisZLOSfOsQSUVFASPVcBohLcGynWlpaThz5gzcbjd0Ol3YJRb1er3wer3weDxwu91wu93weDxgjMHj8Yif8XUVX16v129ZHQqFotILEGJCeIHFqKgoKJVKKJVKcZ1SqRSX/HPf71wscZyccLlcOHfuXKXkeBxZCJjCwkL06tULGzZsEItjaTQavP/++7jvvvtgs9kCZgLUaDQBo5MDuTrdbvclk0iJaHh8HxClGkLynUZdxQNIg1Nx9hN5YCKTYNlOnU6HZs2aoaCgALm5uQE7Wr4uUAfM1/t26r4dNlA5sysXD0C5AOEigr/n4sP3vdfrFUWK2+0WP+dtUalUUKlUfm3WaDSVxELFtvq2n8Pf83b6Ln1fAESBxNvo21beRn5+vqKKb6dQKKBSqfyua0Uh5Lv0vb6BBFWgax6o7fzF28/bG6j9vM0ulwupqal+RTN9CalZ/OSTT3D8+HGxBPrUqVNx+eWX49tvv0XLli39Gnn55ZdjwYIF8Hg8SElJ8avoySksLAyozDp27IgJEybA4/HA4XDAZDLB6/WSgCFqTHo6sGeP8F6qh8T8/PL3UgmYivlnSMBEJsG0ndHR0WjSpEnAjrbii693uVx+2wTyalTs6Dm+nW3FDpm/5501f887eC5S+N8N7cHg+26oY/Drw8UCFxO+19bpdFb6n1R1rX33G+hcAr18BSD/H6hUKrFqtq+o1Wg01V6LkJrF4uJi2Gw2qFQqxMfHI+FChbpA6YWLiorEk2nZsiX279/v9zljDDt27AiYxW/y5MmYPHlypfUkYIiakppa/l6qOkS+Ve2lEjC+GYA1GgrijVQawnb6eiWqinEgggu/5ppLZCw4pHfNQw89FHD99u3b0bVrVzFVtdfrxaJFizBixAgolUqMHz8et912G86dO4eMjAwAwIoVK1BQUFCpDkR1qNXqKotCEYQvF+qvAZBHIUU5CBiphtII6SHbScgRWcje66+/Hu3bt8ezzz4Lm82Gt956C9u3b8f69esBABMnTkTLli0xcOBAPP7448jNzcVLL72EESNGoE+fPjU+TnR0NOxS5WQnwgrfkUmpOm45CBjfh24qIxa5kO0k5IgsHMLLly+HzWbDwIEDcfXVVyMnJwc//PADhg4dCgCIiYnBH3/8gaFDh+LRRx/Fu+++i0cffRRLly6t1VihVquFQ6qc7ERY0aRJ+XupBIxvMUepPL6+8fAkYCIXsp2EHJGFB6ZXr17YvHkzcnJywBhDRkZGJWFiNBrx0Ucf4aOPPqrzcTQaDZxSlfUlwgrfWUhSiQffKcxSiSjf2DxK3RG5kO0k5IgsBAwgBBfx+JaGQqfTwWazNegxiEsD32BVKdImeL3+WXDlEH9yIeaeiEDIdhJyRBZDSKGC/wgpmR1xMaSuOGE2CyKGI1UyPV/CLOcYEUTIdhJyJOIEDAAKRiMuitSFyyuWEZCDB4YETORCtpOQIxElYGIvVOkrkUORGULWSG2nK4YbSBWH4zuURgImciHbSciRiBIwvByB1Xd6B0EEwHf6sBTemIoTPqTywPhO305MlKYNhPSQ7STkSEQJGF44jILRiIvhOwNIivxdFQWMVN4PX8+P79RyIrIg20nIkYgSMDzTL/0IiYtx/nz5e9+U/qGi4hCSVEG8vp4f3+zERGRBtpOQIyRgCCIAvoUUpbhdKgoYqXKw+CavMxqlaQMhPWQ7CTkSUQJGr9cDAEp9xwcIIgC+t4jFEvrjV+wnqqgm3+D41oGSqiYUIT1kOwk5ElECJu7CYyxF0hMXw2Qqfy/FQ2fFY0olHnxjbyiRXeRCtpOQIxElYOgpgqgpOTnl76UQMBUDh6WKgbkwexYAkJQkTRsI6SHbSciRiBIwfCog/QiJi5GXV/5eihp2vtO4AemGkHyF04U+jIhAyHYSciSiBExCQgKUSiXyfHsnggiAr6e8opgIBRUztkuVg8V3GjUNIUUuZDsJORJRAkalUiE5OZl+hMRF8R02krr8i1Ip3Swk3+sgVTZgQnrIdhJyJKIEDCC4QikQjaiOipWgpSY52T8jbiihWbMEh2wnITciTsDo9XoaxyWqpaTEvxq1UoJfie8xpUzhL0X8DyFPyHYSciMiBUyZFKlVibCh4u2hUIS+Db4CJjk59Mfn0E+F4JDtJOSGRI5p6YiNjSU3KFEtvpWoNRppYj/69BHaYbNJU0ySQx4YgkO2k5AbESdg4uPjcebMGambUSu8Xi9KS0tRVFSEnJwc5ObmIiEhAYMHDwYAFBYW4tFHH0VRURFKSkrgcDjgdDrhdDrhcrnwwQcfYNCgQQCAF154Ac8//zy8ASJT09PTce7cOfHvzMxM5ObmQq1WQ6fTITExEfHx8TAYDLjhhhtw/fXXAwBOnTqFNWvWIDo6Gnq9HgkJCYiNjUVCQgKSk5ORmJgIpRTjMHXEagVmzwbeeIOBMSdKS0tx5kwZSkpKcP78eZhMJlx77bXi9u+++y5OnjwJm80Gu90Oq9WKkpISlJWVYcKECZgzZw4AYMeOHRg9ejRcLhfcbrf48nq98Hq92LlzJ3r16gWvF3j88ZlYufJTxMTEIDo6Wry2Wq0WXbp0wWeffSYe/5FHHkF0dDQSEhIQFxeHlJQUNGrUCCkpKWjSpIlYiK8upKcDubkNL2RcLhfUajUA4X5etWoVSktLxVdZmXD9S0pK8PTTT6NDhw4AgJdffhkfffQRXC4X7HY7HA4HXC4XnE4nWrVqhaNHj4rHSE1NRWlpKTQaDWJiYhAbGyu+ZsyYgSlTpgAAjhw5gh9++AFxcXGIj49HUlISEhISYDAYEBsbi/T0dERFRTXsBQky586dw/79+3H27FmYTCYUFhYiLy8PFosFGo0Gn3/+ubjtsGHDcOzYMdhsNvFePXnyZFjaTuLSJuIETFxcHIovRGgOHDgQxcXFiI2NRXx8PGJiYnDXXXdhxIgRAICDBw/izz//RFJSEvR6PWJiYqDVaqHRaKDT6dCqVSsoLowvuN1uKBQKKJVKeL1eeDwe0ZBqNBoxEZTZbMauXbtgs9lEg1xaWgqLxQKLxYK5c+ci6ULGsKlTp+KHH35ASUkJWIXH8JEjR4oChjGGRYsWVXnOvjMHFApFQPECoNL6c+fO4bxvVUMfOnfuLL7fs2cPZs2aVeXxz549i4yMDADAPffcgx07diAxMRGpqakwGo3Q6/WIi4tD165dMXr0aABCh7Z//36xs4mJiYFKpYJarYZarUZ0dLTYifBr4/F44PF44Ha74XQ6YbfbwRgTj+3xePD111+LQs9qtaKwsBAFBQUoLS3FPffcg/Hjx8NqBfbvXwSDYSZcVZSi9ng8oihbvHgxdu7cGXC7Nm3aiO+joqKQ71tkqQKLFzPcfTewfz/Qq5dHFKHFFSKKeUcPAE6nE2+88UaV+/zvf/+LGTNmAAB++uknfPDBBzAYDIiPj4fRaERSUhKSk5ORlJSEsWPHAhDKKJw9ewY6HbB0qXCttVotXK4oKBQKREVFife91+uFzWaDzWZDWVkZrFarKDhsNhuGDRsGzQUX1vfff4+///4bhYWFohA/d+4ccnNzcc0112DJkiUAAJPJhGnTplV5TrfeeqsoYKxWK06cOBFwu4r3M29jWVkZzGYzcnyyFY4fP158/9dff+HRRx+t8vh5eXlISUkBANx11104dOgQkpOTkZKSAoPBIN7PvXv3xpAhQwAI/6dTp05BrVZDo9GIdkStVkOlUvldU243HA4HHA4HbDYbTCYTrFYrdDod+vTpAwBwOBx45ZVXxM8sFgtKSkpgsVhQXFyM1157DWPGjAEALFq0CE8//XTA84n1zVYI4PTp08jOzvZbV1xc7Gc7CUIORJyASUxMhNlsBiAYgH379vl9zo04AGzYsAGzZ88OuB+FQgGPT6Rnjx49cODAgYDbPvjgg3jrrbcAAHv37sXw4cOrbN8dd9whChi1Wg3LhUI83LPRqFEjpKeno3fv3n7ntHDhQiQlJSEuLk40jvzVrFkzcdsnn3wSTzzxhGgwFQoFGGOiF8CXPXv2iE+0ZWVlYsdvsVjQpUsXcbtmzZphxowZsNvtKC0tRXFxMUpKSlBUVASTySSKNwA4dOgQ/vrrr4DnfsMNN4gC5syZM+jVq1eV12nDhg0YOnQoAMED8eabbwbcrkOHDjh48CAAQKlU4vbbb69SlHTrNgrZ2cDnnwNRUWrRK6DT6aDX6xEbG4uUlBSkpqbC7XaLHfPs2bNRUFCA6OhoxMTEiE/qMTExaNKkibj/Ll26ICcnR+y0+EupVEKpVKJvXyV27RK2veeehVi37j1RHNjtdpSVlcFut4uF9QChk37zzTdFAVxcXIzc3Fzk5eUhNzcXyT4BNDt37sTKlSsDnntcXByefbYYK1cCv/8ONGo0GNnZWQG3feKJJ/DKK68AAH7++WeMGjUq4HYAkJ2dLV6Dzz77DCtWrAi4na/IzsjIwNSpU8XrrtfrodPpEBcXh9jYWD/xPGfOHEyfPh0qlUoUWlwkVPSSFBYWit5Ju90u3stWq9VPaLZt2xZz5syBxWKByWSC2WyG2WwWvZu+13TPnj1V3s933HGHKGAOHTqE7t27V3mdduzYgczMTADA3XffjY8//jjgdn379sWff/4JQLifn3/++Sr36est6dy5M4YNG4bGjRvDaDSKtoR7VH1Zu3YtoqKiEBMTI96rBoPBz3YShByIOAFjMBhQVlYGr9eLH374Abm5uaLxt9lsfsKgffv2mDp1Kkwmk/hUyQ2gSqUSn5gA+L0HhKftQIY0OTkZgwcPhk6nE93X/IktLi4ORp+Sv/Pnz8cbb7yB2NhYqKqZRxsVFYU777yzRucfaD/8qbqiwU9LS6vRPnv27In//ve/Ndr2o48+Ql5eHkwmE3Jzc8VrW1EUOZ1OdOvWDU6nU+zAXS6X+PL1QnAPDD8PlUoFjUaD6Oho8UmZf86vU1xcHAwGAxISEsQn5x9+6AyuV8eNmwSb7Xqo1dG42GhBdd4CX9RqNRo1aoQ//gBeeQXYsEFIkuf1At26AePGAa+/DvTsCSQmKnHvvdH48stovPNOIm69NfA+o6Oj8dBDD9Xo+DfeeCO6dOmC0tJSmEwmFBUVobCw8ELHHoWHHxa2UypxQaTZ/IZmPB4PGGN+w4FqtVr0kHHxxgVHTEyMn+fwuuuuQ/fu3WE0GpGWlob09HRRkPt2ojqdzm+IrDqSk5P9BEV1cEF/Mfr27Yu+ffvWaJ+LFi1Cfn6+eC9zD5TZbPazJS6XC61atRLvX7vdLg7xcu8thw8ZarVaaLVa6HQ6xMfHIzY2Fp06dRK3U6vVeO6558R7mdsQPvTVtGlTcduJEydi4sSJNTqn1q1bB1zvazvDaUiYuHRRsIpjE5c48+bNw5NPPonS0lLopCowQ0iC2SwMzxw9KsR1WK3CuoICIU2+QiHkW9m4EcjKApo0EZbBmoVUXAw88ACweLHwt8EATJsG3HUX0LVr5e2vvx747jvg66+BG28MThs4ubnA2rXC+RuNwrDRypXA1KnAVVeVT93mMZs6HS4q5GqL1wu89x5w/jzw7LOAVhvc/RPBhWwnITci0gMDgH6EEQRjwJIlwP33Cx02AHTsCEyZAgwfDlx2GeAzKoNPPgGmTwcGDw6eeNm9WxAkWVlCR/3EE8BDD1Vf44i3NZh5YBgDFiwAHnsMcDqFDL8zZghBy3fdVXn7u+8GvvhCEF1VeYHqwunTQGamIB4B4I47AJ9RHEKGkO0k5EbECRg+RJOfn+83vEBcmjgcQgd9IT4UgwcLT/tDhlQtTnj2WV9RUx82bwZGjxY8Pj17Ch6Vdu0u/j0+nbsek4j8cLkEYcYnnNxxB/Daa4IHJlRtAICtW4Hx48vFC0DiJRwg20nIjYgVMCaTSeKWEA2N2QxMmAD89pswRLRgAXD77Rf3qgRTwOzfD4wcKSSEmzxZ8O7UVAwEUzy43cDNNwPffitci0WLgBtuuPj3+BBShYkqdWbbNmDECGG/LVoAJ08KS0L+kO0k5EaNBczff/+N5cuXY8OGDcjOzobVaoXRaERGRgbGjx+PG2+8EY0aNWrItgYF7ga1Wq0St4RoSAoKgCuvBP7+G2jcGPjxR6CaSSB+8ElKPnHCdcJuF0RLWRlw0018dlPNvx9M8fDYY4J4iYsDfv4Z6NevZt/jiVeDMWJw9qzgeSkpEa7LRx8JQ2inTgnHoVEJeUO2k5AbFw0lP3PmDEaOHIkePXrg008/RfPmzTFt2jQ88cQTuOmmmxAXF4e5c+eicePGmDmz6rwZcoHnPKCMkpcuZWXAmDGCeGnXDvjjj5qLF0CIDQHqn4H32WcFD0ybNsB//1v7INjCQmFZ3RBPTfj+e+CttwRBtmpVzcULEDwvkMcDTJoE5OcLwvKzzwTB0ry5EJdz6lT99k80PGQ7CblRrQfmr7/+wogRI9C2bVts3rwZ/fv3rzRdGBASRK1YsQL33XcfioqKsHTp0oDbyQGeY6XAdwCeuGRgTAhG3b5dGJrYuBG4kMeuxgSj0/77b2D+fGFK8uefCzOOagNjQFGR8P7CLVsnTCYhEBcQ2jNwYO2+Hywx99ZbgpBs3Bj48svy6tqtWwPHjwOHDwuB1YR8IdtJyI1qPTDz5s3DmDFj8Pvvv2PAgAFVipKYmBjcdNNNWL9+PTZu3CjrG5wHn1WXEZUIX/73P2HWjF4veBtqK16A4AybPPaYME343nuFWU61xWoVRIxOV97Z14Xnnwfy8oBBg4AqcjJWi9stLOvThqwsgCeBXbgQ8I3/5Kl/9u+v+/6J0EC2k5Ab1Zqlb775BkVFRdUmUfOlW7duyMvLk3WSI41GA4PBgCL+eEtcMuTnA/fdJ7x/6y3AJ2Frraiv12HTJiHOJC4OeO65uu0jGCIqOxv44AMhaHnBAv8K1zWFJ2euz0/6qaeE2WBTpgizsXwhARM+kO0k5Ea1ZmnHjh1ISUnBrbfeih07dlSqxxNwhzIWLxyDwUCBaJcgjzwiBO8OGyZMna4rvEJEXRO3vfiisHz44boP//DiifUZxnr1VUGMTZpUuxggX+rrgdm3T5g2rtUCL79c+fO2bYXlyZN12z8RWsh2EnKiWrXRu3dvvPvuu9iwYQP69OmD9u3b48UXX8TJMLc2Go0GTv6YTVwS7NghBIZqNMIwRTBCsOqyj+xsIe5GqxUS59UVPpW7rgLGbBamSgPlwzd1ob5ibv58YXnnnYBPSS4RPsTnUwSdkDFkOwk5Ua2AUSqVmD17Nk6dOoVffvkFAwYMwOuvv46WLVtiyJAhWLRoUVhGpEdHR8POIzWJS4KnnhKWDzwgBIbWB+5E9KnVWWOWLBFiV8aPBxIS6t6G+g5jffyxMAw1fHj5ME1d4E7XujhWc3KAr74SvltVuSYuYHJy6na9idBCtpOQEzUyS1FRURg+fDg++eQTnD9/HkuWLIFWq8X06dPRqFEj3HXXXWHllaEf4aXFnj1CzInBIKTory9cNNTlQZNnua1hfccqqc/QDa8xBAiCLhjUpWLahx8KOXUmTgRatQq8TXS0IGLcbqG8ACFvyHYScqLWz1U6nQ5TpkzBZ599hmeffRZOpxP//e9/8fjjjzdE+xoEcoNeWvDOevr0+k055vBhm9ra6aws4OBBITnbVVfVrw318cBs3izElDRrJuTDqQ986Ki23hGvt7xo5cVmP3GP2YkTtTsGEXrIdhJyosbPdy6XC1u3bsXatWuxdu1a7Nq1CwqFAkOHDsVNN92ESZMmNWQ7g4pKpYKbP+ISYU1BgZBXBBCmLAcDXhWZB9LWlDVrhOVVV9U/i299PDC87tOUKfWbPQSUf5/PRqopmzcLyemaNRPqT1VHy5bC9sePCwHYhHwh20nIiWrNo8vlwrvvvotff/0VGzduRElJCdLS0jBixAg8/PDDuOqqq5CcnByqtgaNqKgoeGjA/ZJg0SIh4HXUqPIZLfWF10DigbQ1ZdkyYVlxqnB9qG0gsc0GLF0qvJ86tf7H50Kstgm2P/tMWN5yy8VFFC9seeRI7Y5BhB6ynYScqFbA5OfnY8GCBWjVqhUef/xxjBo1Cj169AiLqdLVERUVBW9tHykJ2cFYeczJzJnB229cnLC0WGr+ndJS4NdfBcExfnzw2lJbAbNundDuXr2Ck9m2LgLG4wFWrBDeT5ly8e15MUeKgZE/ZDsJOVGtgMnIyAir4Fwisti7V8gzYjQKFZ+DRXy8sCwurvl3/vpL6OR79ap/7aL6sHy5sLzmmuDsry4BzX/8IdRxatOmZiIqPV1Ynj9f+/YRBBG51Ck9FWMMRUVFOH78OE6ePIns7Gy0bdsW44P56NmAeL3eGmcXJuTLt98Ky+uvL49bCQZ1ETBbtgjLyy8PXjuA2sWeuN3ADz8I74MlYOrigeHel/Hja+ZBSksTlrm5tWsbEXrIdhJyolZ34q+//orHH38cBw8ehMXHv56amooZM2aEjYDxeDzQBrPHIySBx5xcf31w98tnMtUmY/qGDcJyyJDgtKEus3+2bBE8H23bAp06BacdPB6opjOyGCv/v0ycWLPvpKYKy7y8WjWNkACynYScqJWAmTp1KmJjY/H666+jY8eOaNq0KdLT04N6Q//yyy9wu90YGWBM4Ny5c/j000/hcDhw4403olMdrbTH40FUXVOLErLg5EngwAEgNvbis1xqS2KisDSZara91ysMIQFC0cRgwB9yazPhY/VqYTlhQnAyEQPlVbRrmj3+6FFhOrnRCFxxRc2+YzQKnh6TSUi+V5/6T0TDQraTkBO1isbVarUYPHgw7rrrLgwcOBAtWrQIqnhZtmwZRo4ciQMHDlT6bMmSJWjbti0++OADfP311+jatSteffXVOh3H6/WGfSBypPPTT8Ly6qvrP2W5ItwDU1hYs+0PHRI6+CZNyr0J9aUuU7nXrROWI0YEpw2+7aipB4a34cora15+QKkUrh0AnDlTu/YRoYVsJyEnauWBeeGFF3Drrbdi7NixyMjIwPnz53HmzBnk5OQgJiYGDz/8MBR1fPTbuHEjbrrpJsyZMwdz5szx++z48eOYPn067r33Xrz66qtQqVT4z3/+g3vuuQc33HADWlWV5rMKXC4X1MHu9YiQwnOuBLOz5vASADWdhfT778IymPEver2wLC2t2fbnzgkZiWNigP79g9eO2npgfvlFWF55Ze2Ok5EhJLI7d658WjUhP8h2EnKiVgLmu+++g9vtxhif9J6JiYlITk5GkyZN8NBDD9XJvejxeDBr1izMmjULr776aiUR9NlnnyEtLQ2vv/66uP+77roLr7zyCj755BO89NJLtToe/QjDG4ejvKMcNSr4+4+OFrwCDocQvHqxW2XrVmEZTOFQWwHDBd2wYeVxK8GAD6eZzRff1m4v/79cfXXtjtO0qbCkqdTyhmwnISdqLGAYY/jpp58wZcoUPPPMM4iNjUVycjI0da0258OKFStgsVgwceJE/Pvf/0a7du0wdOhQMdp9/fr1GD9+vJ84UiqV6N69Ow4dOlRpfwcPHsShQ4egVCqh1WqRkJAAnU6HZs2aISEhAW63m36EYcxffwmxEl26lA89BBOFQhgKOn9eePHOtSoaYgZSbKywrGmtVD59ety44LWhtu3YtEnw1HTvHrjydHXwa0xDSPKGbCchJ2osYBQKBf7v//4Pjz/+OG666SaMHTs2aI147bXXkJubi2HDhkGpVMLr9aJz58749ddfYTQace7cOT+vDychIQFnAli8pUuX4rnnnqu0vkWLFjhx4gRsNhuiecEbIuz4809hWdMg0bqQkSGIl5yc6gWM2SwErmq1QM+ewTt+TIzg+XE4BM9GdberxyMk0QPqX/uoIrVJ6vfbb8KyLsN6jRoJS8oFI2/IdhJyosbRWCaTCUqlEv3798fEiRMxe/ZsvP766/j4448xf/587N+//6L7yM/Px8GDB7Fz50788ccfKC4uhsfjwY4dO3DrrbciNzcXTqcTGzZswMmTJ7Fw4UIAQHx8PKwBBuFLSkpg4IP0PlSV6jrmgm/dZrOJ74nwY/t2Ydm3b8Mdg1fIuFgg7+7dwrJbt+AGEysU5d6Pi8Wf7NkjCIxWrYLvkeIzgmoylPXzz8Jy6NDaH4cETHhAtpOQEzX2wBw5cgQLFixAUVERmjVrhk2bNsFsNsNkMiExMRGtW7dGly5dqt3HpEmTsH//fkRFRcFgMODVV1/F8OHD4fF4cMsttyD1whSOoUOHYuzYsdi8eTOefPJJpKWlBfS0HD16FBMmTKi0vmPHjpgwYQI8Hg8cDgdMJhOsVisSLkRnOp3OoAx9EdLApyz36dNwx+CziS7WofK29O4d/DbUtCo2j38ZPjz4beAxMBfLiVNYCOzaJXii6jKtnQsYSmYnb8h2EnKixgKmX79+OHbsWL0OtnHjxkrrGGPQarU4XSF6z263Q38hknHo0KF488034Xa7xbiYrKwsHDhwIGAA7+TJkzF58uSAbWCMobS0NKDnhpA/JpMwWyU6Oji1fqqCd6j5+dVv15BiiguYixWV3LRJWF51VfDbwNP8X0xY/PabkMTu8svrFkRMHhj5Q7aTkBvVDiE98sgjWLBgQY13lpubixEjRqCoFilMFQoFhg0bhsWLF8N1IV/58ePH8fPPP4txNjfddBOKi4tx3333wWq14sSJE7jhhhvQqFEjXF3L6Q42mw0ejwex3D9PhBW7dgnL7t2Dn//Fl5oGr/J4nIYYzqpJLhinU6g9BAQviZ4vfCjtYkJu8+b6tYHqIckfsp2E3KhWwPTp0wcPPvggZs2aBfNF5lEePXoUo0aNQl5eHuJ45F8Nef7557Fz505069YN1113HXr37o1+/fqJXpTmzZvj66+/xvfff4+4uDi0atUKxcXF+Pbbb6GrZdpOXgKhtm0k5MHffwvLXr0a9jjcRlcXvJqTA2RnC9t26BD8NtQkjf+2bcKMrE6dymsKBROeE+didaF4AG9dBUxCgpB92GKpXfI+InSQ7STkRrVDSDfeeCNSU1Nx6623YvHixRg3bhwGDhyINm3aQK/Xo6SkBH///TdWrFiB7du3o0+fPli5cmWti3316dMHWVlZ+OSTT1BYWIjbb78do0aN8ps2PXHiRAwePBhbtmyBVqvF0KFD6zSdjwsxHg9DhBcHDwrLhhw+Aso9D9XV5+HeoMzMmmedrQ38Fq3Ooclz0AS7nAKHPx+UlVW9jdMpVAZXKOruiVIohJICublCPE1GRt32QzQcZDsJuXFRpTF06FAcPnwYX331FZYvX44nnngCZT7WzGAwYMKECXj22Wdx1VVX1TlHQGpqKp544olqt0lMTMS4eia6KL7wKBnPSw4TYcWRI8KyITwevvByAtV5HvgMpB49GqYNKSnCsqCg6m34ENZllzVMG3hCPatVqPkUKIv8gQNCwr+2bcs9V3WBBIy8IdtJyI0auUp0Oh2mT5+O6dOngzEGk8mEkpISGI1G6PX6OpcPkALuBqUfYXhy+LCwbOh083z2TXXTqPlwVkMLmKriQrze8vwvwcwC7ItKJYiY0lLhFUig8Gnt9Z2JVZcq4EToINtJyI3ajfVACLpNSkpCErc2YUbphYQWfIYTET4UFwtP6NHRF8+OW18aNxaWZ89Wvc3evcKye/eGaQPPZltVev0DB4TOvmlTIQdMQ5GQIIiXoqLAAoYXcBw4sH7HqemUbUIayHYSciPiyooWXnikTuTWkggb+Cz+Nm0CD2UEE67Pq4pdLykRMvCq1Q0Xj8Nz0VTlBeKzj/r3F2JIGgoeHBwoHoix8hlItS3gWBH+kzSZ6rcfomEg20nIjYgTMHkXrHBaQ0zZIBqUrCxh2bp1wx/LN3g1UGLnPXuEzrtr1/LpzsHmYrOQGjr+hcOHsgIJmKwsYX1KiiAs64PRKCwvlv2YkAaynYTciDgBYzabodVqKR12GOLrgWlooqKEOkCMBc4Fw2cgBbP+UUUuVpGaez4asiYUUP1Uah7/ctll9fcC8ZlfNIQkT8h2EnKjWgFz+PBh9OvXD1999RXcbneo2tSgWCwWymMQppw4ISwbMt7DFx6rGGhIg8e/NFQAL1DukQiUBbe0VBB0anXDtgEo9wQFygjMS6AFow0UAyNvyHYScqNaAZOamork5GTcfPPNaNKkCR544AHs2LEDjLFQtS/oFBQUhG0AcqRz8qSwbN48NMfjAiJQh7pvn7Ds1q3hjs+nEgcSMP/8IyzbtWvYjMQAwDPHB/IEcQHTqVP9j8MF48WS5hHSQLaTkBvVCpjExET8+OOP2L17N2655RZ8++236NOnD7p06YJXX30VZ6uboiFTioqKYOQ9ExFW8OnEPO18Q1NVOQGns3wKdUMKmOqmFfMEdpmZDXd8TlVDWb4BvMEopcAFzEWSfhMSQbaTkBs1ioHp0aMH5s+fj9OnT+Pnn39Gr1698Pzzz6Np06YYMWIEli1bBq/X29BtDQqlpaU0DTBMOXdOWIYqyRkvuut0+q//5x9hXdu25fEhDYFeL8y2stmERHG+8BichqjBVJGqsvGeOiUMr6WkBGdYj4J45Q3ZTkJu1CqINyoqCldddRU+//xz5Obm4oUXXsCGDRtw7bXX4q677mqoNgYVq9VK1VTDEJdL6NiUyvJZMQ0NH5qpKB4aOgMvR6Go2ivBPTD9+jVsG4ByT1BFYRGsBHac6oaqCOkh20nIjVolsnO5XPjzzz+xdu1a/PTTT9i5cye0Wi3Gjx+PGTNmNFQbg0phYSGN44YhvPNMSmqYukOBqErA7NghLEMxfJOYKHg5iovLhZvVWp6DpiGHsDhVeUa4gAmWF4j3jVZrcPZHBBeynYTcuKiA+eeff/Dbb79h3bp1+OWXX2CxWJCamooRI0Zgzpw5GDt2bFhFppvNZvoRhiF8+ChU8S9AeX6XitWRQ1URGwg8A2jfPiH+pGPH8mGuhqQqYbFli7C8/PLgHIcPx1EMjDwh20nIjWoFzKlTp9C5c2eoVCr069cPjz32GEaNGoUePXpA2dCpUBsAl8sFu92O2PpUnCMkIdTxL4BQsgDwTyTndpcPITVkDhgOFzC+8Sc7dwrLUAgooDyY2VfAuFzl1yFYifT4c1BJSdWFIwlpINtJyJFqBUyzZs3wzz//oGnTppfE2CdVUw1fpBAwgYaQDhwQxESrVuWJ1xqSQAG0XMAEK/bkYnBhcaGWHwAhB43TCbRsWR6nU1+USuFYFoswZEYZ6+UD2U5CjlT7jKNQKNCxY8dLQrwAVIwsnOFTqEMpYLgHwHeCHc+/0tABvJxAwyoHDwrLrl1D04ZAgcTHjwvLtm2Deyw+QkH1kOQF2U5CjkSUk9Z+YSwgmo8NEGEDT+YWyjIsPFjYtxbS4cPCsl270LSBB9AWFAhLxsoFTIcOoWlDoFICPCtyy5bBPVZVM54IaSHbScgREjBEWMCfyEM5rBDIA3P0qLAMlYCp6JE4cUIYYklLK69W3dBwB6xvQr9Dh4RlsOtS8Qf8ijlnCGkh20nIkYgSMDSOG77w4YtQ/uu4gPGtnME9MMEeOqmKitmAQxlAzNFohJw0brfwAhoukLiqqeuEtJDtJORIRAkY84VeMKEh06cSDQL3QIQyk3nF6speb7nnoWPH0LSh4lRuXoOpe/fQHB8QrgOfDWW3C4KOC7nOnYN7rKqmrhPSQraTkCO1SmQX7lAgWvjC4y9CmXKoYs3SEyeELLFpaaETUhVnIfH4l1AJKE5MjNAGm02YTm02C7ExwR7GonpI8oRsJyFHIkrAcDcoPUWEH3zoQhXCO5YH7/KhpG3bhGWfPqFrQ8VEdnwWVLA9HxeDJ8xzOMq9Lx07VvZS1RcqJyBPyHYSciSihpBKLgQSUDKm8IN3aKF8AOQChs9G4iUEQlF/iOObTM/rFfKvAKELIub4FrbkmYgbYhgrUMAwIT1kOwk5ElECxmKxQKlUQsf98kTYIMUspIoChleADmUArW/sydmzwjI1NbRDaUC558vtLo8D6tQp+MehadTyhGwnIUciSsAUFRUhISEhLMsgRDJ2uzCEolKF1gPjO2zl8ZR7YEI5hMSFg8sFZGUJ71u3Dt3xOXyoiLFyL1BDzMSqOOuKkAdkOwk5ElF3Y1lZGT1BhCFFRcLSaAx+zEV18Km8arWQ/6WsDGjWLHT5V/ixeVuys4X3zZuH7vgc3ynlPBdOsHPAAOVBvL5J8wjpIdtJyJGIEjAulwtq3iMQYQMvIhjq4XceOBsTUx642hDDJtXBY2AcDiAnR3gfyorcFXE6haEspbJhhBSvL0VDSPKCbCchR0jAELJHigBeoFw4GQzAkSPC+1AlsOP45kXh9aAaNQptG4DybMQ5OYIXJj293DsUTPj0dBIw8oJsJyFHIkrAuN1uqEI5D5cICtwTEuos5r4ChsefNMSwSXXwAGK3Wz4CBgAaN26Y4/AhJN/K14T0kO0k5EhECRh6ighPeFZW7o0IFb6eH159uVWr0LaBCxivF8jLE96HMgaHw5P6nTkjLJs2bZjj8FlXVAtJXpDtJORIRAkYp9MJDU9oQYQNXEjwHCGhgneiOl35EFKo86/4IqWA4QHNvA1NmjTMcfg0eRpCkhdkOwk5ElEChtyg4QkfTgh17hM+dBUVBZw+LQSutmwZ2jb4ljMoKBCWKSmhbQNQ7gXjM8IaKpA4IUHwwlitNJVaTpDtJORIRAkYj8eDKO6TJ8IG3pGF0gPDWHk9HpNJ+LtZs4YJXL1YOwBh+nh+vvBeCgHDc+Jwz0hDxeEoFOX75jE/hPSQ7STkSEQJGMYYJWIKQ3gAaSgfAK1WwQMTE1Pu+Qh1/AtQng3YahWGcQyG0AczA5UFTFpawx0rI0NYnj3bcMcgagfZTkKORNwdqQhlJjQiKDidwjKU3o9z54Rlejpw8qTwXooEclw48DggPktHqnbwkg48X0tDQB4YeUK2k5AbESdgmG9QAREW+E5nDhW5ucKyUaPyKdShjn8BysUbDyiWqhgwD+Llw2oNKWD4DKfTpxvuGETtIdtJyA0SMITs4Z03n2IbCrgHxlfASDGExIUDF3G82GEoYazcA8MDqhsyDocLRX7dCXlAtpOQGxElYKKiouDhQQVE2MA7z1DO4uTxF02alA8htWgRuuNzKuZDCXU2YqDcC8Tfa7XC1PKGgher5Ll3COkh20nIkYgSMCqVin6EYQj/l4VyCJ7P+ElLKxczDZW8rTr4DCyexE+KXGLcC8RJSmrY/wX3dJEHRj6Q7STkSEQJGI1GAwdPaEGEDb5VoUMFD1ZNSChP3taQM2+qgntg+ASQUA6jcXgAMaehh7G4p+v06XLxSkgL2U5CjkSUgImJiYGNZycjwgYpSgnwWA/GhE40Li70pQwA/2R6QGiH0TjFxf5/N3RCwZgYYfaX211euoCQFrKdhByJKAGj1+tRWvFxkpA9fBglNjZ0x+QZZ7n3R4rkcUC5eOAzsKSYyVqx3wqFF4jHwRw71vDHIi4O2U5CjoSVgGGMISsrCyd5VGUt0el09BQRhkgxjbpiLZ6GnDZcHVzAcNEghReoYmXoULSB15ziNagIaSHbScgRyQXMd999h1atWqFdu3bo0qULMjMzccUVV2DAgAEYMmQI9u3bBwDYu3cvBg4ciNatW6Nly5YYOXIkTpw4UatjqdVqOH2nVBBhgW9V6FCRkyMsedyN0Ri6Y/vC+wwu3qSICak4hBQKDwwXMAcPNvyxiItDtpOQI5JX5+rbty9mz54Nm80Gl8sFl8uFkpISvPfee2jdujVatGiB3NxcDB8+HO3bt8fvv/8Ot9uNRx55BNdccw127dpV4xTXGo2GfoRhCPfAhGoIibHyWUhcwPAqyaHGbheW/Nz5lPJQUnEqdygyynfvLiz//rvhj0VcHLKdhByRXMA0bdoUc+bM8Vt3//33Q6fTYeXKlYiNjcW8efOgVquxdu1a6C88hi9ZsgTt27fHunXrMGLEiBodi/8IGWOUFjuM4EMYoRIwJSVC7IteXy4gpMqAy4/PhZQU9fS4B0avF7xhofjpdOkiLPfvFwQl/VylhWwnIUckFzAV+eOPP/Duu+9i4cKF6NChAwDgp59+wtSpU0XxAgDt2rVD69atsXXr1koC5uDBgzh06BCUSiW0Wi0SEhLQuXNnaLVaMMbgdruhliKhBlEnePr6UHlBePyL0Vh+bKkETMU8MFLMQuJlFeLjBQETimKSjRsL07WLioQ8PE2aNPwxiaoh20nIEcljYHxhjGHu3Lno2bMnpk+fLq7bt28fOnXqVGn7tLQ0nAkwz3Lp0qW49tprMXHiRIwaNQqXX345NmzYgNgLj/CWilGJhKzxzckSCnxjbriAaOipw1XBjx/KAOaK8J8YL7IYiiEkhQLo1k14v2dPwx+PqB6ynYQckZWA2bp1K3799Vc89dRTfnEtKpUqoNvS4/FAparsRAqUMbKwsBDGC5GYJt4jErKnrEwYRtFoGjZ9vS++MTdSzIDyhXuDUlOFZcWsuKGAJ/LjmYhD1YaePYXl7t2hOR5RNWQ7CTkSsiEkr9eLUaNGIT8/HzabDU6nE4sWLcKgQYPEbT744AO0adMGEydOFNcpFAqkpqYin0dV+pCTk4MJEyZUWt+xY0dMmDABHo8HDocDJpMJCoUCiRfGIIp4kg9C9hQUCMuUlNDFQfh6PUIdf1OR8+eFJc9DI4X3nsfA8JlYoRIwffoIy23bQnM8omrIdhJyJGQCRqlU4plnnkFpaSlUKhViY2PRkz9iASgoKMA333yD+fPnI6pCpGLXrl3xxx9/4OGHHxbXnTx5EqdPn0a/fv0qHWvy5MmYPHlypfWbN28GABRXnBdKyBauW0OZh4UPIRkMoZ8B5YvTKYiHqKhy8eD1hr4d3AOTkSEsQyVgBgwQlr//Lpx3KIauiMDEx8cDINtJyIuQBvEO4BYpAJ9++im8Xi9uvvnmSp9df/31uPvuu/HPP/+gU6dOcLlcmDt3LuLj49G3b98aH58HAVNGyfCBd558CCUU8MDd+Pjy+A8pBAwPnk1JKZ8+LYUHhreD14IK1Uyopk2BZs2Emkj795fHxBChh2wnIUdk80zz8ccfY9y4cUgKUCluypQpuOyyy5CZmYmbbroJffv2xddff4358+fDUIvgBHqKCD945xlKAcNvj4SE8vdSBPHyYycllQuYUE+jZqxc0HEvUCiT6fER5k2bQndMojJkOwk5IgsB43a7odfr8cQTTwT8XK1W45dffsHbb78Nt9uN3r17Y8+ePZgxY0atjsMD0Qp4YAUhe6SoBM0fMnU6aYeQeAI5na48H0yoq1E7ncKQkVpdfg1Cmc9s8GBhSQJGWsh2EnJEFnlgVCoVduzYUe02UVFRmDlzJmbOnFnn48THxyM6Oho5PE88IXvOnROW6emhO6Zv9WveWUuRf4VP+IiPl64dvoU0ufcnlB6Y4cOF5fr1wjWQ4v9AkO0k5IksPDChQqFQID09Hef51A5C9vAg3lAOIfEOOiqq/L0UAaR8wofR6C+qpGhDYqIwnASE9lq0bAl07SrMBtu4MXTHJfwh20nIkYgSMIAwHdDMB/UJ2ePbiYcKXwEjRSFJDvfAJCVJ1w4e8nAhBEISeFaF77+Xrg0E2U5CfkScgImLi6NAtDCCx8CEchq1r4Dh1aBDHXsClOegiYsrb0co0vj74nv+PA9PqKdyX3utsFy+XJpiloQA2U5CbkSkgCnhA/uE7OFD7jwHSSjgHbRS6f8+1PAg3piY8nYESDwdkjbo9eXXINQCpnt3oH17QcyuXBnaYxPlkO0k5EbECRij0Yg8/lhPyBqvV5oYGF/RIkXgKifQFG4ehxIq+KSTpCR/z1QoUSiAu+8W3r/7bmiPTZRDtpOQGxEnYBo1aoS8vDywUPcERK3JyxOGDJKSQh+8CggdJ++0Q+35APzz0UgFj8MxGsuFXagFDABMmyZ4gTZuBP75J/THJ8h2EvIj4gRMWloaPB4PCnmVPEK2ZGcLy2bNpDk+Y9JOo5YygJjjWwuKDyFJ4Y2KjwemThXev/pq6I9PkO0k5EdEChgAAYtDEvKCp/Fv3Di0x+XeFre7fMhGCq+Db/I6qYay+BCS0ShdDAzn0UeF/82SJcDhw9K0IZIh20nIjYgTMLz0gJWnWCVky4kTwrJly9Ae11fA8E5bitkvvh4YqTxB/GcSF1eei0aqZHKtWgG33y4IqKeekqYNkQzZTkJuRJyAibsQEWnhvnFCtpw6JSxDPYTkK2D4tGU+nTiU8AkfBkO5NybUsUC+ZRWkLCjJeeYZoS3ffQf89JN07YhEyHYScoMEDCFbTp8Wls2bh/a43MPgdPqLmVDjW0rA5fJvW6jwrcck5ZRyTtOmwPPPC+/vuac8RodoeMh2EnIj4gSMTqcDQGXhw4GzZ4VlkyahPS73cjgc0qXwB8o9MFIKGH5cX68LT2gnFQ8+CPTqBZw8Cdx2m3QxOZEG2U5CbkScgOFPEZSQSf7wGJimTUN7XF8PjJTTqH2z4EqViZd7nlQqabMS+6JSAf/7nyDsli8XhpWIhodsJyE3Ik7AxMbGAqAfodwxmYQZMHq9dB4Y3+rHPIg2VDBWHn9iMEgnHnyDh6USUYFo0wb46ithdtbLLwPz5kndoksfsp2E3Ig4ARNzoQco44P7hCzxnYEU6iELnnfFavUXM6HE6RSGRlQq4eVbViDU7QDkJ2AAYNQo4NNPhfvjySeBWbNC/3+KJMh2EnIj4gSMUqlEdHQ0jePKHCmT2CUmCkuTqTz/Co8FCRUVPS6+OWFCiRyHkHy55RYhL4xWC/znP0CfPsD27VK36tKEbCchNyJOwABCMJpNinmxRI05flxYhnoGEiBknQXKg2iloGIWXh5MLGUMjJQBzdVx883Ab78JeWL27gX69QMmTQL++osCfIMN2U5CTkSkgDEYDJSMSeZwAdOuXeiPzYsnFhdLl0CuoliQagjJt4Cjb0yO3OjbF9i3D3jsMeF/9c03wrpmzYCZM4FFi4AdO8rrSxF1g2wnISckmFshPXq9PuJ+hB4PkJMD5OYChYVC/ozSUsBsFoZKioqEmA+HQ+i0+ausTNjOZhP+drsrP9Xyqs0qlTDdtuJLpxMKMsbFCZ1fQoLgWUhIEGaSREcLL71e+Ds2Fjh2TNh3qMsIAOVDSGazdLOQfIUTY+Udb6jFAy+loFRK14aaotMJdZLuuw94/XXg+++FchQLFwovTmKicI/xe1GnK4+zUiiE/7XLVT6N3mYTroNSWV4fy+sV7g23u/w3wUtPKBTCtjx+KSpKeKnVggDVaIT7XasV2hAdLbQnLq68PYmJQgX2tDRhfePG8vB8RaLtJORLRAoYtVoNVy2DGsxmoHNnwdjExAgGhhvCuDjB6BgMwrr4eGGd0Sj8rdcLn2u1gvGKialdh+h2C+LCZisXExaLMMRRWip0LCUlwjbnzwtVnHNyBGFSXCxsa7GUd0bhhNEY+mNyLwePw9HppBu6OXoUeOst4X9rNJZ7h0LNp58Cq1YJ71u1kqYNNaVJE+Cdd4C33xaGkTZuBHbvBg4cEDx7JlN5ksBwgtsao1EQNGlpguBJSwOSkwU7YzQK28XHC3/HxpbbnttuA379FVi8GBgypG5tqIvtJIiGIiIFjEajgbOW0xVsNuDcuWC2ofwpjD+lAYLI8HiEJ0CnU1jyAM76kpYGNGokGLu4uHIvSGJiuYdEo/F/6XTCdvzJ0betHI+n/MXb7PsqLRU8PCUlwstsFtaZTMLfdrvw4h6hkhJBeKWkAN26Befca0NSkjD0wDMBjxoV+vT56enC/8NiAR5+WFj35JOhLyrZsaMQV/LOO8LfvXoBV18d2jbUFYVCGEbq27d8HWNAfr7gWeT3Ig/pYEx4ud3C/1urLfcOKhTlnket1t/rqFKVe1yUSmE77pHhL/6bttnKf9P8nrfZhLZYLMJ9X1Ym/F5yc4WHkeJiIakj//2cPSv8T2oDbxdQP3tSF9tJEA1FRAqYujxFpKYK7mibTXgVFvp7N7gRMpnKDVFBQXlnbbOVD8+UlZUP0dQEpbLci8NfsbHCS68XOrq4OOHvlBRBqKSlCe/5Z3Fx0iRjqyvcWyRF1le1Wnha37gR+PtvYOrU0LfBaAQOHhSEwwcfCILq3ntD347PPwfGjAFeew04cgR4801pSwnUF4VC+C2HG15vua3JzxdETH6+YG/4sLDVKggf7p21WoVlWVn5LDqjEejate7tIA8MIScUjIXjwEL9GDhwIFQqFTZu3CjJ8RkrH193OATjwmMtFIry8XKNpnzcXOr07YR0uFzCPSNVFWhAOL7dLq8p1ETN4TaGe4zqitS2kyB8CaNn8uChVCohpW5TKMpd0wRxMaSs/sxRKEi8hDM8oL6+SG07CcKXMHYG1x2v1wsFuTQIgiBqBdlOQk5EpIDxeDyICnU0JEEQRJhDtpOQExEpYNxuN1ThFNFKEAQhA8h2EnIiIgWMw+GAVg5ZoQiCIMIIsp2EnIhIAWO32xFNEbQEQRC1gmwnISciUsCUlZVBp9NJ3QyCIIiwgmwnISdIwBAEQRA1gmwnISciUsA4nU5opMwKRhAEEYaQ7STkREQKGApEIwiCqD1kOwk5EXECxu12w+VykRuUIAiiFpDtJORGxAmY0tJSAIBer5e4JQRBEOED2U5CbkScgCkqKgIAJCYmStwSgiCI8IFsJyE3IlbAJCcnS9wSgiCI8IFsJyE3Ik7AWCwWAEBcXJzELSEIgggfyHYSciPiBExxcTEAID4+XuKWEARBhA9kOwm5EXECxmQyAaBxXIIgiNpAtpOQGxEnYKxWKwDAYDBI3BKCIIjwgWwnITciTsDYbDYAQExMjMQtIQiCCB/IdhJyQzYC5vz583jppZcwa9Ys/Pe//wVjzO9zt9uNDz/8EFdccQUGDBiAJUuWwOv11vo4xcXFiIqKomRMBEEQtYBsJyE3ZCFg/vnnH3Tq1AkrV66EzWbDM888g8cee0wUMR6PB2PGjMGDDz6IHj16oFOnTpgxYwZmz55d62OVlJQgNjYWCoUi2KdBEARxyUK2k5AbClbR1SEBt956K86cOYP169dDoVBg9erVGDNmDPLy8pCSkoJPP/0Ud911F7Zt24ZevXoBAJYtW4brrrsOx44dQ6tWrWp8rNtuuw2//fYbTpw40VCnQxAEccnhcrkAAGq1WuKWEISASuoGAMCxY8fQsWNHUdlzF6XZbEZKSgq++eYbXHfddaJ4AYDx48fDaDTi+++/xyOPPOK3v4MHD+LQoUNQKpXQarVISEiATqdD165dsXjxYng8ntCdXJBhjKG4uBiFhYUoLi5GaWkpiouLYTKZUFhYiJKSEjgcDjidTjidTrhcLpSVlaG0tBQ2mw1OpxNut7vSNVAoFIiKioJKpYJGo4FarYZKpYJarYZarYZOp0NSUhLi4uIQGxuL+Ph46PV6JCQkID4+HtHR0YiOjoZer0d8fPwla+TcbjfMZjOsVitKS0thsVjEa2uz2WC322G1WlFSUoKysjLx5XQ64XA4YLfb4XK54Ha7xZfX64XX6xU9jvx3wK+777XVarVQq9UwGAyIj49HfHw84uLiEBcXJ75PTU1FfHx82D4pl5SUoKioCKWlpeKrrKwMJSUlKCkpEa8vf8+vqd1uh8PhgMvlgtPp9LvHFQqFeG9rNBrExMQgNjZWfPlev4SEBCQkJIjvExMTL4n72eFw4Ny5czCZTCgqKkJubq54/9rtdvFedTgc4j3N71WPxwOv14tu3brh9ddfl/pUCAKATATMjTfeiAcffBCJiYlo2rQpXnvtNQwZMgStW7cGAGzZsgXz5s3z+05UVBSaN2+OkydPVtrf0qVL8dxzz1Va7/V68cADD2D//v2IiYlBQkICkpKSxA45JiYGBoMBiYmJojFLSkpCUlIS9Ho9VKrgXC6v1wubzYaSkhJYLBaUlZXBYrHAYrHAarUiNzcXubm5OH/+PAoLC8XPTCYTcnJyYLfbq92/QqEQDTU31nq9HjExMdBqtYiKikJUVBQUCgUUCgUYY/B4PHA4HHC73aLw4cXbuAgym801jjuKjo5GQkICjEYjDAYD9Ho9kpKSkJycLHYMqampMBqN0Ov1YgfCO46YmJigd8BOpxP5+fkoKioSO7/CwkIUFhaKHaHVaoXJZILFYkFxcTFKSkrETtRqtaKgoKBWsVcxMTGIiYmBRqOBVqtFdHS0KA75S6lUii9AEKn8HsnNzRWFUVlZmdhZO53Oao+r0WiQmpqKlJQUpKamIj09HWlpaUhLS4NOp0NCQgKSk5ORmJiI5ORkJCQkwGAwiG2oL4wxOBwOUTxzEcLFd05ODs6fPy8uz58/j6KiIvF/URO0Wi0MBgNiYmKgUqkQHR0tCjyNRiPe44AwDG2320Vhb7fbxd8fD06tDp1OB4PBgNjYWPGaGo1GJCUlQafTISUlBcnJyeK9Hh8fj8TERFEMBeO6MsbgdDpRVlYGq9UKi8WC/Px8mEwm8W9+TvyhJicnB/n5+cjLy0N+fn61++fxLVqtVrQXvvdqVFQUysrK6n0eBBEsZCFgpk2bhrfeegtvvPGGaMhff/11KJVKMMZgtVqRkJBQ6XsxMTEBO/NAHpbo6GjRmNntdpjNZhw4cABmsxklJSU18sqo1WpotVpoNBrodDrx6Vir1Yo/cKVSCa/XC4/HI3Y0LpdL7AB5J3QxoqKikJqaitTUVMTGxiI9PR0dO3ZEo0aNkJ6ejuTkZNELEh8fj6SkJCQmJiIuLg4qlapBnr69Xq/4JGw2m1FaWgqz2Yzi4mLY7XbY7XbRI8SfoouKikRvxb59+1BUVASLxQKHw3HR89fr9aIA450U9wgplUpRiPHOwePxwOPxiCKMt8npdMJqtdaoY+SdO/duxMbGIi0tDXq9HrGxseL/RK/Xi+u4secv3tFFR0cHTRBUxOVywWKxwGw2ix1XcXExiouLkZubi7y8POTl5aGgoAA5OTnYv38/8vLyxGGAQCgUClE8chGgVqvFe5wLAqVSCYVCIXqOnE4nbDab2LHyp/eLjU4rlUqkpqYiIyMD6enp6Nq1K5KSkpCRkQGj0QidTideZ51OJ3r/DAYDDAZD0LwiHo/HT7CazWbxuprNZphMJtFOlJSUIC8vD6dOncKOHTtgNpsv2qnz66rX68Xryu0IFwhRUVFiW/g97HA44HA4YLPZRK9fTUb8VSqV+DCQlpaG9u3bo3///mjcuDEaN24sCte0tDTEx8eLdkytVoet146ITEIqYD777DMUFBSIQxm33347WrRogX/9619wOBxYu3YtMjMz8cwzz2DKlCmIj4/HmDFjkJycLCZR8qWoqAhGo7HS+o4dO2LChAmiV8FkMkGr1YIxhrfffrvSj5QxhrKyMthsNvEJvLi4GBaLBQUFBTCZTOITJB+e4e5W7rbmblbGGFQqld+TIHf586dFbpC5+5p7IOLi4mAwGJCSkgKj0Si5MXG5XGjdujXUajXi4uKwe/du8TzS09Prte+ysjLk5eWJ15Z3vr4dstVqFTtH7nngLy4S+TUHID5186ECPvSi0WhgMBiQlJQkPinzjjAxMREpKSnQ6/UNKjiCjVqthtFoDHj/V4XX6xWHDPgwAvdA+V5/PnTAxTe/x/m15i/uPeLeSy7e+P3N73X+N7/PjUajKASluN4ulwstW7ZEbGwsEhISsHXrViQmJtY5QZvX60VBQYHoPfId1jWbzSgqKhIFP79/ubDmnk7u1fO9h7VaLbRarSiqDQYDoqOjRdvBr2VSUhIMBoMo8AJ5Lxlj2LRpE86ePYuTJ09iwoQJ9b6OBCE1IQvi9Xq9uOeee1BaWgqVSoXY2Fg88MADaNSoEZKSkrBo0SLcfPPNAIQfW6dOndCrVy988cUX6NOnDy6//HIsWLBA3J/VakViYiK++OILTJo0qcHb73K5sHPnTnEYoGPHjg1+TKkxmUxISkoCAOj1+hq79sMZi8WC3r17w2w2w+PxoLCwUOomNTgOhwPr1q0Thzu6dOkidZMalEi8ryPxnIlLn5B5YJRKJT788MNK67Ozs+F0OtG4cWNxnUKhQFJSklg8bMyYMfj444/x+uuvQ6vVAgD+97//we1247LLLgtJ+w8ePIjLL78cANChQwccPHgwJMeVkuzsbPF906ZNJWxJ6Dh58iSOHj0KQPg/RwKHDx/GuHHjAETGvR2J93UknjNx6SO5vzwjIwNt2rTBc889h0OHDqGwsBBffvkltm3bJro5p0+fjuLiYowcORKrVq3CvHnzMGvWLEyePBnNmjULSTt9x7kjJZU2nTOd86VIpJ0vEJnnTFz6SB7EGxUVhSVLlmDmzJnisIxOp8Pjjz+OadOmARCeGLZu3YpZs2Zh3Lhx0Ol0eOSRR/DUU0+FrJ1ms1l8HynFzOic6ZwvRSLtfIHIPGfi0kdyAQMA/fr1w+7du8U8Jk2aNKk0w6Bz587YvHkzbDYbtFptyIP/fGfN8GGsSx06ZzrnS5FIO18gMs+ZuPSRhYABhLiX5ORkJCcnV7udVIXEfHN/8CmPlzp0znTOlyKRdr5AZJ4zcekjGwEjdzp06IDnn38eHo8nYoI76ZzpnC9FIu18gcg8Z+LSRxa1kAiCIAiCIGqD5LOQCIIgCIIgagsJGIIgCIIgwg4SMDWAMYYvvvgCAwcOxOWXX44PPvig2poy4YLJZMKUKVPwySef+K13OBx47bXX0K9fPwwdOhQrVqyoVIPlwIEDmDRpEjIzMzFr1iycO3culE2vE6dOncITTzyBSZMm4bnnnqtU3K4m5/3PP//4nffZs2dDeQq1Zv369Zg9ezYeeeQRrFmzptL5lJWV4cUXX0Tfvn1x1VVXBdxm9+7duOaaa5CZmYn7778feXl5oTyFOvP111+jZ8+eOHDggN/6/Px83HfffcjMzMS1116L3bt3V/rumjVrcOWVV6JPnz546aWXUFpaGqpm14oDBw5g8uTJuPXWW3HjjTfiuuuuw6RJkzBx4kS/3C+nTp3CtGnTkJmZialTp4rJGjmMMXz11VcYNGgQLrvsMrz33nsXLRhKEFJDAuYiMMYwdepU3HbbbWjXrh169+6Nxx57DDfddJPUTasXhw8fRr9+/fDll1+iqKhIXG+z2TBgwAC89NJLuOKKK9C0aVNcd911eOmll8RtVq5cie7du+Ps2bMYOXIktmzZgm7duvll+5QbK1euROfOnbFs2TIolUosWLAAmZmZYkp1ft4vvvii33m/+OKLfvvo1q0bzpw543fep0+fluq0qmXu3Lm48sorsXv3bmzfvh2jR4/2y51ksVjQt29fvPnmmxg0aBCSk5MxevRovPPOO+I2X3/9NTIzM2E2mzFixAisWbMGPXr0QG5urhSnVGP279+P2267DXv27MGJEyfE9cePH0fHjh2xatUqjBgxAlarFZmZmfj555/FbZ555hmMGjUKiYmJGDx4MBYsWIChQ4fK8qFFoVDg66+/RlZWFgCIKSb69u0LjUYDANi+fTs6duyIXbt2YeTIkTh69Ci6deuGPXv2ABBs3LRp03DLLbegdevW6Nu3L+bOnYsbbrihRsUjCUIyGFEtP/30EwPA1q1bJ677/fffGQC2bds2CVtWP6ZPn84GDRrEEhMT2csvvyyuf+WVV1hcXBw7duyYuO79999nOp2OFRYWMrvdztLT09kdd9zBPB4PY4wxh8PB2rdvz+67776Qn0dN6devH3vssceYy+VijDG2f/9+BoAtXbqUMRb4vP/973+H9Xn/8MMPbOPGjeLfixYtYjqdjrndbsYYY3PnzmXJycksOztb3GbevHksISGBlZSUMIvFwhITE9n999/PvF4vY4yx0tJS1qxZMzZ37tyQnkttcDgcrGfPnmzs2LEMAFu5cqX42dixY1nPnj2Z1WoV102ZMoVlZmYyr9fL/vnnH6ZQKNjChQvFz0+dOsVUKhX7+uuvQ3oeNeHgwYMMANu1a1fAz71eL+vevTsbOXIkczqdjDHGPB4PGzx4MJs4cSJjjLGff/6ZAWA//fST+L0///yTAWC///57w58EQdQREjAX4Y477mBDhw6ttL5Dhw7s4YcflqBFwcHj8TCv18uSk5PZvHnzxPW9evViDz30kN+2JSUlTKVSsc8//5ytX7+eKZVKdvr0ab9tXnrpJZaWlhaStgeDzZs3MwBs/fr1jDHhvB988EG/bfh5f/bZZ2zDhg1hf95Llixh0dHRooBp27Yte+aZZ/y2yc/PZwDY8uXL2fLly5lGo2H5+fl+2zz++OOsTZs2IWt3bXnmmWdYfHw8O3HihJ+AsVgsTKVSsW+++cZv+w0bNjAA7Pjx4+zFF19kLVu2FEUqZ9SoUezaa68N2TnUlF9//ZUBYF988QWbMmUKGzJkCJszZw7Ly8tjjDF25MiRgA9bn376KVMoFMxut7M777yTDRw4sNK+u3Tpwh544IFQnAZB1AkaQroIW7ZswaBBgyqtb9WqFU6ePBn6BgUJnsnYt0ptaWkpdu3aVel8DQYDUlNTcfLkSfz+++9o2bJlpYJwrVq1Qm5uLmw2W2hOoB5kZ2dj5syZaN26NQYOHIiysjLs3r074HmnpaXh5MmT2LJlS1ied0lJCV5//XWMGjUKt912G5555hlERUUhLy8PR48exeDBg/22T05ORmxsrPi/7tSpU6Xkkq1atcKpU6f8kqPJhe3bt+Pll1/G/PnzkZ6eXukzt9td6f/cqlUrABD/zwMHDqyU6Vuuv/ecnBwAwJQpU1BcXIxOnTph6dKluOqqq+D1erFlyxZotVr07dvX73utWrUCYwzZ2dn4/fffL0kbR1z6kIC5CFarFQkJCZXW63Q62O320DcoiBQXF8Pj8SAtLQ1AecG36s63uusB+KcslyPLli1Djx494HA4sHr1aqjVapSVlYExdkme95EjR7Bw4UKsWbMGTZo0wZgxYwBAjP2p6zm7XC7ZxUdYrVbceuutuPLKKzF9+vSAnwOVz5n/D/k5x8fHV/quXH/vPBB90aJFWLlyJd5//32sXbsWf//9NzZt2iSej0Kh8PtexXO+VG0ccWlDmXgvgtFohMlkqrS+qKgITZo0kaBFwYMH73IBw4u8VXW+RqMRLperys81Go1sK916PB7MmTMHCxYswN13343XXntNbGtCQgIUCkW15+12u8PyvDMzM3H06FEcOXIEN910E66++mpkZWXBaDQCqPy/ZozBZDLBaDTCarVWec7x8fGyS0l/55134vDhw2jbti3uu+8+sX1LlixBTEyM3znzex4o/x0YjcZqf+/8+3LihhtuQKdOnTBs2DBxXadOndCoUSP8888/SE5OhslkAmPMT8TU9JwvVtqFIKSEPDAXoVmzZjh06JDfOq/Xiz179qB3794StSo4cI8L73xVKhUaN25c6XxPnTqFoqIi9O7dG82aNcOZM2fEp1nO7t270aNHD6hU8tTE7733Ht577z0sW7YM//73v/0ER1Xnffr0aRQWFob1eXPatWuH9957D3l5edi1axfi4uIQFxdX6Zz/+ecfOJ1O8ZyzsrIqTafdvXu3LO/95s2bY9y4cXC73Th48CD++usvAMCff/6JH3/8Ec2aNQOASue8e/duqFQqdOvWDc2aNcPhw4cr7Vuu55yWluYnXgBBhJaUlIAxhmbNmsHlclUaCtq9ezcyMjKQnp4e0MYxxmR7zgQhIl34TXjw73//mxkMBmY2m8V1q1atYgDYX3/9JWHL6g8Pcty+fbu47q677mLdu3f3C2J86aWXmEajYVarlRUWFoqBrRyr1coyMjJkHfDXq1cvdvvtt1f5+cyZM6s875KSkrA7b4fDwZ5++mm/+/bvv/9mANiOHTsYY4xNnjyZ9e/fX5xhxJgQoGswGJjT6WRnzpxhANiyZcvEz00mE0tKSmJPP/10yM6lrtjt9kqzkLp06cJmzZol/u31etno0aNZnz59GGOMrVmzhgHwm422d+9eplAoxBlrcsJisbAzZ874rVuxYgUDwPbv389cLhczGo1+Mw1dLhfr0aMHu+666xhjjC1cuJDpdDpmMpnEbfjsy61bt4bkPAiiLpCAuQhFRUWsUaNGLDMzky1btoy9++67TKfTsWHDhvkZ/nBjx44d7F//+hcDwKZOncp+++03xphgrFUqFZswYQJbvXo1mzt3LlMoFGzOnDnid6dPn87i4+PZ22+/zVasWMG6d+/OdDodO3r0qFSnc1ESExPZsGHD2IwZM9j111/PxowZw26//XZ24sQJxljV5+07I2vGjBksLi5OPO8ePXrI9rzLyspYeno6u/LKK9kvv/zCVq5cybp168Z69uwp3rfbtm1jCoWC3XjjjWz16tVszpw5DAB79tlnxf1MmjSJGY1G9v7777Nly5axDh06sPj4eL+p13KlrKyskoD5+OOPmVKpZA8//DBbvXo1u/766xkA9u233zLGGHM6naxbt26sTZs27Msvv2SLFy9mKSkprF27dszhcEh1KlXy6quvMqPRyH788UeWlZXFFi1axOLj49mVV14pbvPcc88xrVbLXnzxRfbjjz+yIUOGMKVSyf744w/GGGNms5llZGSwXr16se+//5699957TK/Xs0GDBoW1jSMufUjA1ICsrCw2atQoBoCp1Wp27733soKCAqmbVS+effZZ1qtXL9azZ0/Ws2dP9s4774ifbdu2jfXo0YMBYPHx8eyll17yM952u5298MILLCYmhgFg/fv3l/2T2t13380yMzPZmDFj2I033shuv/12NmzYML/cF5faee/cuZP17NmTAWAA2IgRI9jJkyf9tvn1119Zp06dGABmNBrZ66+/Lk6zZkwQAXPnzmUajYYBYEOHDq0y54jccLvdrEOHDuzPP/8U13m9Xvbll1+yRo0aMQCsVatW7KuvvvL7Xl5eHps6dSpTKBRMoVCwm2++mZ06dSrUza8RpaWlbNq0aUypVDIATKVSsSlTprCioiJxG7fbzd577z0WHx/PALBu3bqxNWvW+O3nxIkTbMyYMaKNu/vuuytNnycIuUHVqGuBw+GASqWSXfBiQ2Gz2cTMnoHwer1wOByIiYkJccsalkvtvIuLi6ucZcUpKytDTExMpdkqHI/HA5fLhejo6AZqZWhhjMFms1V7zk6nEwqFAmq1OsStqz05OTk4c+YMWrRogZSUlIDb8HPmM5AC4XA4EBUVJfuYLoIAABIwBEEQBEGEHTQLiSAIgiCIsIMEDEEQBEEQYQcJGIIgCIIgwg4SMARBEARBhB0kYAiCIAiCCDtIwBBEGHH+/HlZ748gCCJU0DRqgvCBMYZ58+ZBoVCgUaNGUKlUcLvd8Hq90Ov1SEhIQEpKCnr16oVNmzbBbrdj5MiRIWnbiy++iFdeeQVFRUVBycficrmQkJCARx55BM8//3wQWkgQBBE6SMAQhA9erxeDBg3C9u3b4XK5qtzu3LlzmDBhAkpLS3HgwIEGb9eaNWswZswYfPXVV5g0aVLQ9vvdd99h0qRJWLlyJUaPHh20/RIEQTQ0JGAIIgCMMdjtdvzf//0fFi1ahJMnT6K0tBRmsxkA0KJFC1gsFqjV6gbPyOt2u9G2bVuMHj0a77//fr33xxjzyz57//33Y8WKFTh27FhYZJ0lCIIAKAaGIAKiUCgQExODlJQUlJSUQK1WIyEhAS1atECLFi0AAGazWezwt27disceewwAsHnzZqxatQqMMRw7dgzffPMN8vLyKh3j8OHD+P7777F9+3Z4vd4q27JhwwacPHkSjzzyiLjum2++wYIFC8AYw6pVq7BlyxYAwN9//41vv/0WpaWlfvsoKSnBl19+iUmTJsFgMGD58uXiZ4888ghOnz6NdevW1elaEQRBSAEVvCCIakhOTkZJSQlcLlcl70S3bt3w4Ycf4qabbsKKFSuwdu1aJCUl4cknnwQAjB8/Hj/++CM8Hg+aNm2KgwcPQq/Xw+v14r777sMHH3yAhIQEmEwmDBkyBOvXrw9Yf2nVqlXo0aMHWrZsKa773//+BwDYu3cvPv74YwDAxIkTRWHSr18/bN26FQqFArt27cJVV12F4uJiDBgwAHFxcfjrr78wceJEAECzZs3Qu3dvrFq1ioaRCIIIG8gDQxDVYDQaAUAcOuJ4vV4UFxeLf6tUKpw+fRovvPACli1bhq5du2LTpk348MMP8euvvyI7Oxvbt28HACxbtgyffvoptm7diqKiInz33XfYtGkTzpw5E7ANe/bsQWZmpt86lUqFjRs3Yvny5diyZQv0ej22b9+O5cuXY9GiRfjzzz9x9uxZAMDjjz+O9PR0HD16VGzThAkT/PaXmZmJPXv21OdSEQRBhBQSMARRDYmJiQAAk8nkt76kpAQAEB8fDwBISkqC2WzG4sWLMXHiRKSmpmLatGmYMWMGBg4cCKPRiOPHjwMAFi5ciNmzZyMpKQm33347brjhBsycORNNmzYN2IZz584hPT29UrtKS0uxatUq9O/fH0ajEU8++SQmTJggihN+PI1Gg9zcXDz11FP46quvMHr0aPTt29dvf+np6cjJyanPpSIIgggpJGAIohp4gG5ZWZnfepvN5vc5YwwZGRm44YYbAAjCxmq1AhDiaVq0aIGsrCwAwJEjR/Dzzz+jQ4cOyM3NxebNm/Hhhx/6Bdb6olKpKn3GGMNll12Gyy67TDweF1WJiYmIi4sTj7do0SLce++9yM7Oxi233IKrr746YMyNSkUjygRBhA8kYAiiGvR6PQBUCorlsSpcCBQWForDTQCQlpYmDuEAQJMmTUSPiFKpRHx8PHbu3InVq1fjiiuugMfjwdGjRwO2ITk5Gfn5+X7ranO81NRUPPfcc9i8eTPWrVuHTZs2YcOGDX77KygoQHJycg2uCEEQhDwgAUMQ1cA9LBaLxW89D+h1Op0ABI+MVqsVP09NTfWLacnIyMDp06cBAFOnTsWhQ4dw4MABHD58GOvXr8c111yDKVOmBGxD586dsX79evhmPAh0PF8B07hxY5w+fRp//vknRo8ejR9//BG7d+/G4cOHAaCSB+aXX35B586da3hVCIIgpIcEDEFUQ0xMDC6//HJ0797db31sbCwGDx6M5s2bAwB69OghDuEAQEJCgp/XpFWrVqLAePLJJzF16lTMmjULHTp0wKhRo6BQKPDFF18EbMPkyZNx9OhRHDp0SFx3seO1bNkSZ8+eRZMmTZCXl4exY8eiV69emD17Nm677TZcddVV4rZHjhzBoUOHMHny5LpcIoIgCEmgRHYEEQQ8Hg8sFosY9Ltv3z4sW7YM//rXvwAIQcD79+/HwIEDxe8wxlBQUID4+HhoNJoq9+31etG8eXMMGzYMixcvBgA4HA44nU7ExsYCANavX4+jR49i1qxZAIDs7Gzk5+ejV69e4t+FhYVo2rSp39ATAMyYMQOrV69GdnY2oqKignRFCIIgGhYSMAQRBmzYsAEjR47EPffcg7feeqvKgN/awBjDo48+infeeQdr1qzB8OHDg9BSgiCI0EBDSAQRBgwbNgxr165Fdna2GHdTX9xuN06ePEnihSCIsIQ8MARBEARBhB3kgSEIgiAIIuwgAUMQBEEQRNhBAoYgCIIgiLCDBAxBEARBEGEHCRiCIAiCIMIOEjAEQRAEQYQdJGAIgiAIggg7/h/abslVDdSCfQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def run_LIF(pars, Iinj, stop=False):\n",
    "  \"\"\"\n",
    "  Simulate the LIF dynamics with external input current\n",
    "\n",
    "  Args:\n",
    "    pars       : parameter dictionary\n",
    "    Iinj       : input current [pA]. The injected current here can be a value\n",
    "                 or an array\n",
    "    stop       : boolean. If True, use a current pulse\n",
    "\n",
    "  Returns:\n",
    "    rec_v      : membrane potential\n",
    "    rec_sp     : spike times\n",
    "  \"\"\"\n",
    "\n",
    "  # Set parameters\n",
    "  V_th, V_reset = pars['V_th'], pars['V_reset']\n",
    "  tau_m, g_L = pars['tau_m'], pars['g_L']\n",
    "  V_init, E_L = pars['V_init'], pars['E_L']\n",
    "  dt, range_t = pars['dt'], pars['range_t']\n",
    "  Lt = range_t.size\n",
    "  tref = pars['tref']\n",
    "\n",
    "  # Initialize voltage\n",
    "  v = np.zeros(Lt)\n",
    "  v[0] = V_init\n",
    "\n",
    "  # Set current time course\n",
    "  Iinj = Iinj * np.ones(Lt)\n",
    "\n",
    "  # If current pulse, set beginning and end to 0\n",
    "  if stop:\n",
    "    Iinj[:int(len(Iinj) / 2) - 1000] = 0\n",
    "    Iinj[int(len(Iinj) / 2) + 1000:] = 0\n",
    "\n",
    "  # Loop over time\n",
    "  rec_spikes = []  # record spike times\n",
    "  tr = 0.  # the count for refractory duration\n",
    "\n",
    "  for it in range(Lt - 1):\n",
    "\n",
    "    if tr > 0:  # check if in refractory period\n",
    "      v[it] = V_reset  # set voltage to reset\n",
    "      tr = tr - 1 # reduce running counter of refractory period\n",
    "\n",
    "    elif v[it] >= V_th:  # if voltage over threshold\n",
    "      rec_spikes.append(it)  # record spike event\n",
    "      v[it] = V_reset  # reset voltage\n",
    "      tr = tref / dt  # set refractory time\n",
    "\n",
    "    # Calculate the increment of the membrane potential\n",
    "    dv = (-(v[it] - E_L) + Iinj[it] / g_L) * (dt / tau_m)\n",
    "\n",
    "    # Update the membrane potential\n",
    "    v[it + 1] = v[it] + dv\n",
    "\n",
    "  # Get spike times in ms\n",
    "  rec_spikes = np.array(rec_spikes) * dt\n",
    "\n",
    "  return v, rec_spikes\n",
    "\n",
    "\n",
    "# Get parameters\n",
    "pars = default_pars(T=500)\n",
    "\n",
    "# Simulate LIF model\n",
    "v, sp = run_LIF(pars, Iinj=201, stop=True)\n",
    "\n",
    "# Visualize\n",
    "with plt.xkcd():\n",
    "  plot_volt_trace(pars, v, sp)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Resources:\n",
    "\n",
    "1. Integrate and fire model of neuron - https://www.youtube.com/watch?v=tNyh56m16GY\n",
    "2. Introduction to neurobiology playlist - https://www.youtube.com/playlist?list=PLs6UM3rH5tQwuhURavq5o4CGtAdGYnqZ4 (Specifically Lecture 2.1 to 2.5)\n",
    "3. Equivalent Circuit Model of the Neuron - https://www.youtube.com/watch?v=LdumhvDBpzQ&t=291s\n",
    "4. Resistor Capcitor Neuron Model - https://www.youtube.com/watch?v=KXnHxZdn8NU&t=425s\n"
   ]
  }
 ],
 "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.10.12"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}