examples/example-0.6-custom-model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Minimal model implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook demonstrates how to implement your own model in `neurolib`. There are two main parts of each model: its class that inherits from the `Model` base class and its `timeIntegration()` function. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# change to the root directory of the project\n",
    "import os\n",
    "if os.getcwd().split(\"/\")[-2] == \"neurolib\":\n",
    "    os.chdir('..')\n",
    "\n",
    "%load_ext autoreload\n",
    "%autoreload 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "try:\n",
    "    import matplotlib.pyplot as plt\n",
    "except ImportError:\n",
    "    import sys\n",
    "    !{sys.executable} -m pip install matplotlib\n",
    "    import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model equations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example we will implement a linear model with the following equation:\n",
    "\n",
    "$\\frac{d}{dt} x_i(t) = - \\frac{x_i(t)}{\\tau} + \\sum_{j=0}^{N} K G_{ij} x_j(t)$.\n",
    "\n",
    "Here, we simulate $N$ nodes that are coupled in a network. $x_i$ are the elements of an $N$-dimensional state vector, $\\tau$ is the decay time constant, $G$ is the adjacency matrix and $K$ is the global coupling strength."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We first create a class for the model called `LinearModel` which inherits lots of functionality from the `Model` base class. We define `state_vars` and `default_output` so that `neurolib` knows how to handle the variables of the system. Next, we define `init_vars` in order to use the autochunk integration scheme, so we can save a lot of RAM when we run very long simulations. \n",
    "```\n",
    "class LinearModel(Model):\n",
    "    state_vars = [\"x\"]\n",
    "    default_output = \"x\"\n",
    "    init_vars = [\"x_init\"]\n",
    "```\n",
    "Next we define a simple parameter dictionary called `params`. In here, we can define all the necessary parameters of the model and change their values later. In this example, we set the timescale $\\tau$, the coupling strength $K$, the integration time step `dt` (in ms) and the duration to 100 ms.\n",
    "```\n",
    "params = dict(tau=10, K=1e-2, dt=1e-1, duration=100)\n",
    "```\n",
    "We are now ready to set up the constructor of our model! This method is supposed to set up the model and prepare it for integration. All the magic happens in the background! We pass the `self.timeIntegration` function and the parameter dictionary `self.params` to the base class using `super().__init__()`.\n",
    "```\n",
    "def __init__(self, Cmat=np.zeros((1,1))):\n",
    "        self.params['Cmat'] = Cmat\n",
    "        super().__init__(self.timeIntegration, self.params)\n",
    "```\n",
    "That wasn't too bad, was it? We are finally ready to define the time integration method that prepares all variables and passes it to the last function that will crunch the numbers. Here we prepare the numpy arrays that will hold the simulation results. We have to prepare them before we can execute the numba code. \n",
    "\n",
    "```\n",
    "def timeIntegration(self, p):\n",
    "    N = p['Cmat'].shape[0]\n",
    "    t = np.arange(1, p['duration']/p['dt']) # holds time steps\n",
    "    x = np.ndarray((N, len(t)+1)) # holds variable x\n",
    "```\n",
    "Next, we make use of a neurolib convention to prepare the initial conditions of our model. If you remember, we defined `init_vars` above in order to use the autochunk feature. The autochunk feature will automatically fill this parameter with the last state of the last simulated chunk so the model integration can be continued without having to remember the entire output and state variables of the model indefinitely. In this line, we check whether x_init is set or not (which it will be, when we use chunkwise integration). If it is not set, we simply use random initial conditions using `rand((N, 1))`. Remember that the convention for array dimensions is `array[space, time]`, meaning that we only fill in the first time step with the initial condition. \n",
    "```\n",
    "# either use predefined initial conditions or random ones\n",
    "x[:, :1] = p.get('x_init') if p.get('x_init') is not None else rand((N, 1))\n",
    "```\n",
    "We're ready to call our accelerated integration part and return the results 🚀!\n",
    "```\n",
    "return njit_integrate(x, t, p['tau'], p['K'], N, p['Cmat'], p['dt'])\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numba time integration "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remember to put this function outside of the class definition, so we can use use `numba` acceleration to greatly increase the performance of our code. We first have to let `numba` know which part of the code to precompile. We do this by simply placing the decorator `@numba.njit` in the line above the integration function. Easy way of getting 100x faster code! ❤️ `numba`! \n",
    "```\n",
    "@numba.njit\n",
    "def njit_integrate(x, t, tau, K, N, Cmat, dt):\n",
    "```\n",
    "Next, we do some simple math. We first loop over all time steps. If you have prepared the array `t` as described above, you can simply loop over its length. In the next line, we calculate the coupling term from the model equation above. However, instead of looping over the sum, we use a little trick here and simply compute the dot product between the coupling matrix `G` and the state vector `x`. This results in a `N`-dimensional vector that carries the amount of input each node receives at each time step. Finally, we loop over all nodes so we can finally add up everything.\n",
    "\n",
    "```\n",
    "for i in range(1, 1 + len(t)): # loop over time\n",
    "    inp = Cmat.dot(x[:, i-1]) # input vector\n",
    "    for n in range(N): # loop over nodes\n",
    "```\n",
    "In the next line, we integrate the equation that we have shown above. This integration scheme is called Euler integration and is the most simple way of solving an ODE. The idea is easy and is best expressed as `x_next = x_before + f(x) * dt` where `f(x)` is simply the time derivative $\\frac{d}{dt} x_i(t)$ shown above.\n",
    "```\n",
    "x[n, i] = x[n, i-1] + (- x[n, i-1] / tau + K * inp[n]) * dt # model equations\n",
    "```\n",
    "We're done! The only thing left to do is to return the data so that neurolib can take over from here on. The outputs of this simulation will be available in the `model.outputs` attribute. You can see an example time series below. \n",
    "```\n",
    "return t, x\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numba\n",
    "import numpy as np\n",
    "from numpy.random import random as rand\n",
    "from neurolib.models.model import Model\n",
    "\n",
    "class LinearModel(Model):\n",
    "    state_vars = [\"x\"]\n",
    "    default_output = \"x\"\n",
    "    init_vars = [\"x_init\"]\n",
    "    params = dict(tau=10, K=1e-2, dt=1e-1, duration=100)\n",
    "    def __init__(self, Cmat=np.zeros((1,1))):\n",
    "        self.params['Cmat'] = Cmat\n",
    "        super().__init__(self.timeIntegration, self.params)\n",
    "        \n",
    "    def timeIntegration(self, p):\n",
    "        p['N'] = p['Cmat'].shape[0] # number of nodes\n",
    "        t = np.arange(1, p['duration']/p['dt'] + 1) # holds time steps\n",
    "        x = np.ndarray((p['N'], len(t)+1)) # holds variable x\n",
    "        # either use predefined initial conditions or random ones\n",
    "        x[:, :1] = p['x_init'] if 'x_init' in p else rand((p['N'], 1))\n",
    "        return njit_integrate(x, t, p['tau'], p['K'], p['N'], p['Cmat'], p['dt'])\n",
    "\n",
    "@numba.njit\n",
    "def njit_integrate(x, t, tau, K, N, Cmat, dt):\n",
    "    for i in range(1, 1 + len(t)): # loop over time\n",
    "        inp = Cmat.dot(x[:, i-1]) # input vector\n",
    "        for n in range(N): # loop over nodes\n",
    "            x[n, i] = x[n, i-1] +\\\n",
    "            (- x[n, i-1] / tau + K * inp[n]) * dt # model equations\n",
    "    return t, x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Running the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We prepare a \"mock\" connectivity matrix, simply consisting of 12x12 random numbers, meaning that we will simulate 12 `LinearModel`'s in a network."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "Cmat = rand((12, 12)) # use a random connectivity matrix\n",
    "model = LinearModel(Cmat) # initialize the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "That's it, we are finally ready to run the model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "model.run()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot outputs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Activity $x$')"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEGCAYAAABhMDI9AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8QVMy6AAAACXBIWXMAAAsTAAALEwEAmpwYAABUTklEQVR4nO3dd5xU1f3/8deZXrZXtrCwLCy9I13B3gJGY69JTPJNoqnGxCQmMeab+kuMJkYT9WuCxlhCjGKviKCAgPTet7C9zU4v9/z+mAFZ2F12kd0F9vN8POYxO3fO3Dl3Bva955x7z1Faa4QQQoiDTH1dASGEECcXCQYhhBBtSDAIIYRoQ4JBCCFEGxIMQggh2rD0dQU+raysLD148OC+roYQQpxS1qxZU6+1zm7vuVM+GAYPHszq1av7uhpCCHFKUUrt7+g56UoSQgjRhgSDEEKINiQYhBBCtCHBIIQQog0JBiGEEG1IMAghhGhDgkEIIUQb/TYYQvtaaHltLzLtuBBCtNVvgyFS6aV1SQWGJ9zXVRFCiJNKvw0Ga34SAOEqXx/XRAghTi79Nxjy3ABEqrx9XBMhhDi59NtgMDksmDMcRKTFIIQQbfTbYIB4qyFyQIJBCCEO16+DwZbnJtoQwAjH+roqQghx0ujXwWDNSwINkWppNQghxEH9OxjyEwPQ0p0khBCH9OtgMKfZUQ6LnJkkhBCH6dfBoJSKD0DLmUlCCHFIvw4GAFt+PBi0IVNjCCEESDBgzUtCRwyiDYG+rooQQpwUJBgODkBLd5IQQgASDFhzXGBSEgxCCJHQ74NBWUxYc1xEDsiZSUIIARIMQLw7KXzAK2szCCEEEgwA2AqTMVojsjaDEEIgwQCAtTCxNkOFdCcJIYQEA/HJ9DBBuLK1r6sihBB9ToIBUFYz1ly3tBiEEAIJhkOsBUlEKltlAFoI0e9JMCTYCpMxfFFizaG+rooQQvQpCYYEW4EMQAshBPRyMCilLlJKbVdK7VJK3dXO80VKqcVKqbVKqQ1KqUt6q27WPDeYFREZgBZC9HO9FgxKKTPwF+BiYBRwnVJq1BHF7gae01pPBK4FHuq1+llMWAfIALQQQvRmi2EqsEtrvUdrHQaeAS47oowGUhI/pwIHerF+2AqSCFfIFdBCiP6tN4OhACg/7HFFYtvh7gFuVEpVAK8C32hvR0qpryilViulVtfV1Z2wCloLk9DBKLHG4AnbpxBCnGpOtsHn64B/aK0LgUuAJ5VSR9VRa/2I1nqK1npKdnb2CXtzW0EyIAPQQoj+rTeDoRIYeNjjwsS2w90KPAegtV4OOICsXqkdYM11gUURLpcBaCFE/9WbwbAKGKaUKlZK2YgPLi86okwZcC6AUmok8WA4cX1Fx6AsJmz5SRIMQoh+rdeCQWsdBW4H3gC2Ej/7aLNS6l6l1PxEsTuALyul1gNPA5/XvTwSbCtKIVzZio4avfm2Qghx0rD05ptprV8lPqh8+LafHvbzFmBWb9bpSLZBybBME6nyYRuY3JdVEUKIPnGyDT73OVtR/GzZUJmnj2sihBB9Q4LhCJZUO+ZUG+EyGWcQQvRPEgztsBWlEJYWgxCin5JgaIetKJlYU4hYqyz1KYTof/ptMBjaYFfTrnafOzjOIK0GIUR/1G+D4eH1D3P1y1fTHGw+6jlbfhKYFSEZZxBC9EP9NhjOKzqPiBHhlb2vHPWcsiYudJMWgxCiH+q3wTA8YzgjM0by4q4X233eVpRMpMKLjsmFbkKI/qXfBgPAZUMvY2vjVrY3bj/qOVtRCjpiEKny9UHNhBCi7/TfYChbyaX7N2A1WXlh1wtHPW0bnLjQbZ90Jwkh+pf+GwzVG0hb8Vfm5kzhlT2vEIlF2jxtSbVjTrcT3tfSRxUUQoi+0X+DofQiAD6rUmgKNbGkYslRRezFqYT2emRFNyFEv9J/gyFtIOSOYWblFrKd2e12J9kHp2L4IkTrA71fPyGE6CP9NxgASi/CUraCeYPOZ1nlMuoD9W2ethUnxhn2SneSEKL/6PfBgI7xWVM6MR3jpd0vtXnakuXElGQlvFcGoIUQ/Uf/DoaCyeDKonj/KiZkT+D5nc+3GU9QSiXGGaTFIIToP/p3MJhMUHoh7HqLK4dezj7PPlZVr2pTxDY4hVhziGhzsI8qKYQQvat/BwPEu5OCLVxoTifZlsy/d/y7zdP24lQA6U4SQvQbEgwlZ4PZhmP3O1xWchlvl71NQ6Dh0NPWAW6UwyzdSUKIfkOCwZ4Mg2fDjje4qvQqoka0zamryqSwD0ohJBe6CSH6CQkGgNKLoX4HQ6IGU3KnsHDHQgz9yeR59iGpRGsDsnCPEKJfkGAAGHFp/H7ri1xVehUV3gpWHFhx6Gl7SRoAod3NvV83IYToZRIMAKkFUDAFtizivEHnkW5PbzMIbc1PQjksBHc1910dhRCil0gwHDRqPlStw+ap4rNDP8vi8sXU+muBxDhDSaq0GIQQ/YIEw0Ej58Xvt73MlaVXEtMxFu5YeOhpx9A0Yk0hog0yb5IQ4vQmwXBQxhDIHQtbFlGUUsSZBWfy3PbnDk3HfXCcISitBiHEaU6C4XCj5kP5Smit5oaRN9AQbOCN/W8AYMl2YkqxEZJxBiHEaU6C4XAj5wMatr7EzPyZFKcW89SWp9Bao5TCUZJGaHcL2pD1GYQQpy8JhsNlD4fMYbB1EUoprh9xPZsaNrGhfgMA9qFpGL4IkRp/H1dUCCF6jgTD4ZSKdyft+wB89cwvmU+SNYmntj4FxIMBkO4kIcRpTYLhSKMvBx2DLS/gsrq4fNjlvLXvLWp8NVhS7ViynIR2NfV1LYUQosdIMBwpdwxkj4CN/wHguhHXEdMxntvxHAD2YWmE9rSgI0ZnexFCiFOWBMORlIKxV0LZh9BczsDkgcwZOId/b/83wWgQx/AMdMSQ2VaFEKctCYb2jLkyfr8p3mq4edTNNIWaWLR7EfYhqWBRBHdId5IQ4vTUq8GglLpIKbVdKbVLKXVXB2WuVkptUUptVkr9qzfrd0hGcXzupE3xK5+n5E5hdOZoFmxegLaAfUgawe2NfVI1IYToab0WDEopM/AX4GJgFHCdUmrUEWWGAT8EZmmtRwPf7qn6LNlRx3eeXddmjec2xl4F1RuhdhtKKb4w5guUtZaxuHwxjtJ0onUBoo2y3KcQ4vTTmy2GqcAurfUerXUYeAa47IgyXwb+orVuAtBa1/ZUZWo9Qf67tpJ15c3tFxh9OSjToVbDeUXnUZhUyOObHsdemgZAcIe0GoQQp5/eDIYCoPywxxWJbYcrBUqVUh8opVYopS5qb0dKqa8opVYrpVbX1dUdV2UuHDMAm8XEovUH2i+QnAvFc2DjQtAas8nMLaNvYWP9RjbEtmDOcBDcLuMMQojTz8k2+GwBhgFzgeuAR5VSaUcW0lo/orWeorWekp2dfVxvlOKwcvbwbF7eUEWsoykuxl4FTXuhYjUAlw29jHR7Ov/YsgBHaTqh3c3oqJy2KoQ4vfRmMFQCAw97XJjYdrgKYJHWOqK13gvsIB4UPWLe+HzqWkOs3NPQfoFR88HqgnX/BMBpcXLdyOtYUrGEpoIgOiynrQohTj+9GQyrgGFKqWKllA24Flh0RJkXiLcWUEplEe9a2tNTFTp3RC5um7nj7iR7Moy6DDY9D+H4/EjXDr8Wh9nBAt9z8dNWt8k4gxDi9NJrwaC1jgK3A28AW4HntNablVL3KqXmJ4q9ATQopbYAi4E7tdYd/Dn/6TltZs4flctrm6oJd9QlNOEGCHlg28sApDvSubL0Sl4sW4Qe5CCwtbHjM5uEEOIU1KtjDFrrV7XWpVrrEq31LxPbfqq1XpT4WWutv6u1HqW1Hqu1fqan6zR/Qj4tgQjv7+hgEHvQLEgbBGv/eWjTF8Z8AbMys8z9MbHGIFGZbVUIcRo52Qafe93sodmkuawddyeZTPFWw973obkMgBxXDpcPu5xHgvGwCGzusUaNEEL0un4fDDaLiYvH5PHWlhp8oWj7hSZcF79f9/ShTbeOuZVmayu16R4CWyUYhBCnj34fDACXTywgEInx+qbq9gukFUHxWbDuKTDiYxF5SXlcVnIZr1mXEKnwEm0J9WKNhRCi53Q5GJRSDyilVE9Wpq+cMTidQZku/r2mvONCE2+E5v2wb+mhTV8a+yWWJ68HILhFWg1CiNNDd1oMrcAipZQbQCl1oVLqg56pVu9SSnHlpEJW7GmkrKGDgeSR88GZDqsfP7SpMLmQCSPOoNJWi6ej1oYQQpxiuhwMWuu7gaeB9xKB8F2g3RlST0VXTC5EKfjPxxXtF7A64oPQ216G1k9C4Mvjv8yK5A1E9rRiBDsYoxBCiFNId7qSziU+yZ0PyAK+qbVe2vmrTh0FaU5mlWSxcE0FRkdTZEz5IhhR+PjJQ5sGpQzCOjINszZRtXZ3L9VWCCF6Tne6kn4M/ERrPRe4EnhWKXVOj9Sqj1w1pZDK5gAr9nYwXpBZAkPOhjX/ACN2aPPlc6+l3tLM3hWbeqeiQgjRg7rTlXSO1npZ4ueNxNdV+N+eqlhfuGDUAJLtFhau7qA7CeCMW8FTATveOLQpLymPusF+CmrT2VsjrQYhxKntuE9X1VpXAeeewLr0OafNzGfG5/Pqpipag5H2C5VeDMn5sPr/2myeMGcmVm3lnXde6oWaCiFEz/lU1zForQMnqiIni2vOGEgwYvDC2iMnfk0wW2DyLbDrHWj8ZH6/rKGF+F1h0vfY2Na4rZdqK4QQJ55c4HaE8YWpjM5P4amVZR1PjjfpFjCZYeUjhzYppUibUMBk30j+tvLhXqqtEEKceN05K+kbSqn0nqzMyUApxY3TB7GtupU1+ztYoS0lD8Z8DtY+CcFP1mNIm1iARVtgh4+Pqj7qpRoLIcSJ1Z0WQy6wSin1nFLqotP1KmiA+ePzSbJb+OeK/R0Xmv51CHvh4ycObbIWJmFKt3G+bya/X/17DC2ruwkhTj3dvcBtGPB/wOeBnUqpXymlSnqobn3GbbdwxaQCXt1YTYO3gzmQ8ifEp+Re+TeIxS9sU0rhGpfD2NahHKgt56XdMhAthDj1dGuMQcc73asTtyiQDixUSv2uB+rWp26cPohwzGDhmk5OXZ1xG7SUw7ZPAsA9MQelFdca8/jT2j/hj8haDUKIU0t3xhi+pZRaA/wO+AAYq7X+GjAZ+FwP1a/PlOYmM3VwBv/6qKzjK6FLL4L0Ylj+l0ObrAPcWPPcXOqbQ62/lgVbFvRSjYUQ4sToToshA7hCa32h1vrfWusIgNbaAD7TI7XrYzdML2J/g58lHa3uZjLD9K9BxSoo/2Sw2TUpB2u1wdWZl/P3TX+nzt/B64UQ4iTUnWBwaK3bjMYqpX4LoLXeekJrdZK4eEweuSl2/m/Z3o4LTbghPuvq0vsObXKNzwEFnzeuJGJEeHDdg71QWyGEODG6Ewznt7Pt4hNVkZORzWLi5hmDWbarnq1VnvYL2ZNg2tdgx2tQHZ8ryZxiwz4sHcuWIDcMv4H/7vwvm+s392LNhRDi+B0zGJRSX1NKbQSGK6U2HHbbC2zo+Sr2rRumFeG0mnm8s1bDtK+ALQmWfdJqcE/KIdYc4tbUG8h0ZvK/K/5XTl8VQpwSutJi+BcwD1iUuD94m6y1vrEH63ZSSHPZ+NzkAl5cd4C61g5OXXWmxyfX2/xfaIhPoucYlYmymdEbWrljyh1satjE8zuf78WaCyHE8TlmMGitW7TW+7TW12mt9x92a+yNCp4MvjirmHDM4MlOL3i7DUxWWPZHAEw2M85xWQQ21nFx3oVMzp3M/R/fT3OwuXcqLYQQx6krXUnLEvetSinPkfc9X8W+NyQ7iXNH5PDUiv0EI7H2CyXnwqSbYf0z0BK/9sE9dQA6bBDYUM+Pp/0Yb9jLn9b+qRdrLoQQ3deVFsPsxI/3ACO11sla65SD9z1au5PIrWcW0+ALd37B26xvxu8TZyjZBiZjHeDG91E1w9KHcf3I61m4YyGb6mVBHyHEyas7ZyUlAW8qpZYqpW5XSuX2VKVORjOGZDJhYBp/XbKbSKyDQeS0onir4eMF0LQPpRTuqQOIVHoJV7Ty9fFfPzQQHTM6aHkIIUQf685cST/XWo8GbgPygCVKqbd7rGYnGaUUt589lIqmAIvWHei44FnfA5MFlsRnCXFNzAGLCd+qapJsSfzgjB+wuWEzT219qpdqLoQQ3XM86zHUEp8rqQHIObHVObmdOzKHEQOSeei9XR1Pk5GSD2d8CdY/DfU7MTktuMZl4V9XhxGKceHgC5lTOIcH1z1IRWsn3VJCCNFHujNX0teVUu8B7wCZwJe11uN6qmInI6UUXz97KLvrfLy+ubrjgrO+DRYnLP4VkBiEDsUIrK9DKcXd0+/GpEzcu/zejhcDEkKIPtKdFsNA4Nta69Fa63u01lt6qlIns0vH5lGc5eYvi3d1/Es9KTs+h9Lm56F6E7ZBKVgHuPB+eACtNQPcA/j2pG+zvGo5i3Yv6t0DEEKIY+jOGMMPtdbrerAupwSzSfG1OSVsPuDh3W21HReceTvYU+HdX6CUImlmAZFqH+G98RXfrh5+NZNyJvG7Vb+jPlDfS7UXQohjO57rGDz97TqGI10+qYCiDBd/eHNHx2MNznQ48zuw43XYswTXxGxMLgveD+ID1yZl4mczf0YgGuDXK3/di7UXQojOdfk6hsOuX0g5Xa5j8NQf33TYVrOJb583jC1VHl7b1MlYw7SvQWoRvPljlDk+1hDY0kC0MQjAkNQhfG3813hz/5u8tve146qLEEKcaN0ZfP5tV7adKlb+9zkWfO82PPWddAd14rIJBQzLSeK+t7YT66jVYHXA+fdA9UZY/zTu6XmgwLui6lCRL4z5AuOyx/GLFb+gxldzXHURQogTqVen3VZKXaSU2q6U2qWUuquTcp9TSmml1JTu7L87hs88C20YvPHwA2ij+7Oemk2K755fyu46Hy+srey44OgroPAMeOcXWJxRnKOz8H1UjRGOX+BmMVn41exfETWi/PTDn8pZSkKIPtdr024rpczAX4iHySjgOqXUqHbKJQPfAlZ2dd/HIy13AHNv+RJlm9az9o1XjmsfF40ZwJiCFO5/ZwfhaAfhohRc+CvwVsOHfyZpVj46GMW/+pPWwaCUQdwx+Q4+PPAhz2x/5rjqIoQQJ0pvTrs9Fdiltd6jtQ4DzwCXtVPuF8BvgWA39n1cxp5zIcUTp7D0qb/TUFne7dcrpbjjguGUNwZ4+qOyjgsOnBpvOXzwAPY0D7ZBKbS+X4E+bGqNq4dfzayCWdy3+j72tnSy9oMQQvSw3px2uwA4/LdvRWLbIUqpScBArXWnf8Irpb6ilFqtlFpdV3f86ykrpbjgf76JxW7n9b/chxHr/vxFc0uzmVmSyf1v76DFH+m44Pk/j9+/fhfJcwuJNYfwr/+k7kop7p15L3aLnR8u/SGRWCf7EkKIHtSdwecFSqm0wx6nK6UeP1EVUUqZgPuAO45VVmv9iNZ6itZ6SnZ29qd636T0DM770m1U797Jiuef7fbrlVL8+NKRNAciPLh4Z8cF04pgzvdh28s4TKux5LpoXVKBPmzgOseVwz0z7mFzw2Ye+PiB4zkcIYT41Loz+DxOa9188IHWugmY2I3XVxK/evqgwsS2g5KBMcB7Sql9wHRgUU8OQB80fMZsRs6ey4r/PEPFlu5PiT06P5WrJhfyjw/3sb/B13HBGbdDVinqtTtJnp1LtMZPcFvbhtd5g87jmuHXsGDLApaUL+l2XYQQ4tPqTjCYlFLpBx8opTIASzdevwoYppQqVkrZgGuJj1sAh7qssrTWg7XWg4EVwHyt9epuvMdxO+9LXydtwABe+dPv8Htauv36Oy4YjtVs4jevbeu4kMUGl/4BmvfjankCc5qd1vfKjzoT6c4z7mR4+nB+/MGPqfZ1cp2EEEL0gO4Ewx+A5UqpXyil/hf4EPhdV1+stY4CtwNvAFuB57TWm5VS9yql5nen0j3B5nTxmW/fRcDbyut/ua/bp7Dmpjj46pwSXttUzUd7Oxl+KT4Lxl6FWv5HkifZCZe1EtrTNojsZju/n/N7wrEwP3j/B0SN6PEckhBCHJfuzJX0BHAFUANUAV8h3t3TZVrrV7XWpVrrEq31LxPbfqq1PmomOa313N5qLRyUM3gIc2/6EnvXrWH1y//t9uu/fOYQ8lId3LNoM9GOFvMBuOCXYHHiLv8JpmQrnrf2H9VqGJw6mJ9M/wkf137MQ+se6nZdhBDieHV3PQY7MAi4E/g58b/8TyvjL7iE0mmzWPr0Aiq3dW8CWafNzN2XjmJLlYcnV+zvuGByLlz4S1T5ElKGlBPe5yG0q/moYvNK5nH50Mt5dOOjLC5b3M0jEUKI49OVC9xKlVI/U0ptA/4MlAFKa3221vrBHq9hL1NKccFXv0lqdi4v/fHXtDZ2b+bTS8YO4KzSbP7w5g5qPJ1cijHxRig5F/fuOzAnm2l58+hWA8CPp/+YUZmj+OGyH8r1DUKIXtGVFsM24BzgM1rr2VrrPwOn9YLFdpeby+68m3AgwKI//IpoONzl1yqluHf+aMIxg1+83EmLQymY9wDKZJDsfoVIeetRZyhBfLzh/rn3YzPZ+Pbib+OLdHLWkxBCnABdCYYriI8pLFZKPaqUOhdQPVutvpc1cBAX3/Zdqnft4J3HH+7WHEaDs9zcNncoL2+oYunOTi7ASxsIF9yLu+lhzO5wfKyhnQn58pLy+P2c37Pfs5+7l90t8ykJIXpUV658fkFrfS0wAlgMfBvIUUo9rJS6oIfr16eGTZvJ9M9dy6bFb7Huze7Np/TVuUMoznLz0xc3E4x00sCa9HnUkFmkxB4lcsBHYFP7XVdT86byncnf4e2yt3l046PdqosQQnRHd85K8mmt/6W1nkf84rS1wA96rGYniZlXXs+QyVNZ/I9H2Lf+4y6/zm4x88vPjmFvvY/73trRcUGTCT77MC77Cqy2alpe24vuYEK+m0fdzCXFl/DntX/mzX1vdvdQhBCiS7p7VhIQv+o5MS3FuSe6QicbZTJxye3fI6uwiJf++Gvq9nd9AHjm0CxumFbEo0v3sGZ/U8cFUwtR8+8nVT9IrCmE98MD7ddFKe6ddS/js8fzo2U/YmPdxu4ejhBCHNNxBUN/Y3e5uPyue7A5XTz/m3u6dabSDy8ZSX6qkzsXru+8S2n0Z3FMGYPDtArP23uJ+dqfRM9utvPA2Q+Q5cziG+9+gypvVbvlhBDieEkwdFFyZhaX/+BnhAN+/vubnxMO+Lv0uiS7hd9+bhx76o7RpQRw0W9JzXwLHY7heb3jspnOTP5y7l8IxULc9u5teMPe7hyKEEJ0SoKhG3IGD2Het++ivnw/L/3xN8SiXZsae/awLK6fVsRjS/ewZn8n02XYk7Be+0vcljfxraonUt3xL/yStBL+MPcP7Gnew/eWfE+m6RZCnDASDN00eMJkzv/y7exb/zGvPngfhtG1Szp+dMlICtKdfOuZdXiCnfwSL5hEykXDUPhpfvK9Tk9NnZk/k5/N+BkfHPiAuz+4G0N3f4lSIYQ4kgTDcRh7zgWcdeMX2bF8KW8/+pcuXVeQZLfwwLUTqWoJ8qPnN3b6GvOsz5M6aDOhhlQCb77X6X4vH3Y535r0LV7d+yq/W/U7ucZBCPGp9dtg8DQE2LLsAEY7F5R1xRnzrmD6Fdew8d03WfLk/3XpF/KkonS+e34pL2+o4t9rKjouqBTuL3wVq62c5ve8GNWdzLsE3DrmVm4adRNPbX2KxzY+1t1DEUKINvptMGz7sIrF/9zGv3+9isodnZxK2omZV9/IxIvmseaVF1i+8Okuvearc0qYPiSDexZtZnddx2MIypFE+tUTMXQKLY8/C+GOB7uVUnxvyveYN2Qef1r7J/6949/dPhYhhDio3wbDGZ8p5oIvjSboi/DCfWt5/W8b8dQHurUPpRRn3/JlRs85j+UL/8WH//7XMV9jNinuv2YiNouJ2/+1lkC44zEK25iRuEdG8XmmEn7qbuhkjQiTMvHzWT/nzIIz+cXyX/DS7pe6dSxCCHFQvw0GpRTDpuRywz3TmTa/mP2bG3jqnhUs/+9uwsGuL4yjTCYu+Oo3DoXDB88+ecxupQGpDv549QS2VXv48X87H29IvWYOJkeMph0T0O/+utP9Wk1W7pt7H1MHTOXuD+7mlT3dm8ZDCCGgHwfDQRabmSmXFHPDz2cwbEouH7+xn3/+dAVbPuj6+IPJZObCr36TsedcwIrnn2Xp0wuOGQ5nj8jhW+cO4/m1lfyzk7UbTA4L6ddMIKKL8Sw+ABue63S/DouDP5/7ZybnTuZHy37E63tf79IxCCHEQf0+GA5KSrdz3udHceVdU0jNcrL4yfj4w4F2FtBpjzKZOP/LtzPuvItY9eJCljz52DHD4ZvnDOOcETnc+/KWTqfMcI7MxDUxi9bY1YSf/wPsX97pfp0WJw+e8yATsidw19K7eGv/W106BiGEAAmGo+QOTuGKOydx/q2jCHoj/Pf3H/PmY5tobexk0Z0EZTJx3pduY8KFn2HNKy/yxsP3E4t23C1lMin+ePUE8lKdfP2pNZ0u7JM2fximZBuN0e+i/3UD1HS+upzL6uKh8x5iXPY4vr/k+9JyEEJ0mQRDO5RSlJ4xgOt/Pp0zLh3MnvX1/OtnK1j1yl4inQwWH3ztOV/4H2ZceT2bl7zDoj/8kkio41/4qS4rj9w8GW8wyq0LVuEPtx8kJqeF9M8NJxrNpyV8HfzzCmgu67Qubqubh85NhMP73+f5nc8f++CFEP2eBEMnrDYzU+cN4fp7pjF4XBYfvbSXf92zgp2razrtJlJKMfOq6znvS19nz9rVLPzfnxDwtnZYfsSAFB68fhJbDnj45tPriHUwtuEckYF72gC8gYsI+kvgycvB1/mEfkm2JP56/l+ZWTCTn334M57Y/ETXDl4I0W9JMHRBSqaTC788hsvvmIjDbeXNxzbzwn1rqSvv+Jc9wPjzL2Het39AzZ6dPPuzH+Cpq+2w7NkjcvjZvNG8vbWGX7+6tcNyqZcOwZLrotG4i1iTB566EoKeTuvhtDj589l/5vxB5/P/Vv8/Hl7XvRXphBD9iwRDN+QPS+eqH57B3BuG01jl47lfrWLxU9sItHa8JnTp9Nlc8cN78TY28NSPv8uBHR3/0r9l5mA+P3Mwjy3by5PL97VbxmQzk3ndCIyIicb0R9BVm+Cfn4NQ5yFlNVv53Vm/47KSy3ho/UP8+qNfE+viPE9CiP5Fnep/OU6ZMkWvXr2619835I+w6pV9bFxcgcVuZupnihkztwCzuf2sbags54Xf3ktrYz0X/s83GXnm2e2Wixmarzyxmne31/Knaycyb3x+u+W8Kw7Q/MJuUicHSN56HRRMgRsXgj2503ob2uAPq//AE1ue4OyBZ/Pbs36L0+Ls3sELIU55Sqk1Wusp7T4nwfDpNFX7WPbvnZRtbiR9gItZVw1j0OjMdssGWj0suu9XVGzZxPQrrmHmVTegTEcHSSAc45bHP2JteROP3jyFucNzjiqjtabxX9sIbKon61wPjg9ugoFT4YaFYE86Zr2f2voUv/3ot4zJGsOfz/kzmc726yyEOD1JMPQwrTX7NzWw7N87aakNMGhsJrOvHEZaruuosrFohLcfe5hNi9+kZMo0Lvr6d3C4j/5F7glGuPZvK9hT7+Wft05jyuCMo8oYoRi1D63DaA2Tc0ENlje+GA+H658FR+ox6/1O2Tvc9f5dZDmzeOi8hyhOLT6+D0AIccqRYOglsajBhsUVrHplL7GIwbhzBjLlksHYnZY25bTWrH39JZY8+X+kZOUw/44fkT3o6F/K9d4QV/91OXXeEE9/eTpjCo7+ZR+tD1Dz4Dos6Xay5+7H9OKXIWcE3PhfSMo+Zp031G3gG+9+g6gR5fdzfs+M/BnH/wEIIU4ZEgzt8IejaA1uu+XYhbu7b0+YFS/uZuuHVTiTrEz/bAkjZ+ShTKpNucptW3j5/t8Q9Pk470tfZ/Scc4/aV2VzgKv/uhxvKMo/b53G2MKjwyGwrZGGBZtxjssmY3IF6rkbIbUAbnoB0gYes77lnnK+8e432OvZyx2T7+CmUTehlDrm64QQpy4JhnYs+HAfv319GxePyePKyYVMK87AZDqxvwxr93tY9txOqna3kF2UzJlXDyNvaFqbMr7mJl554HeUb9nIuPMuYu4tX8Zqs7cpU97o59pHVtAajPDkrdMYP7DtPgA8i8vxvLGP5HMGkjqiGp66Oj7WcNN/IXv4Mevqi/j40dIf8W75u8wvmc9Ppv8Eh8XxaQ5fCHESk2Box6bKFv65Yj8vb6jCG4pSmO7kc5MK+dykQooyjx4bOF5aa3aurmH587vxNoUYdkYuM68oISn9k1+6RizGsmefZNWLC8ksLOLSb955VNdSRZOf6x5dQbMvwhO3TmViUfpR79P8/C58q6pJv2IY7qI6ePIKiIXh2qdg8Oxj1tXQBn9b/zceWv8QozNHc//Z9zPAPeDEfBBCiJOKBEMnAuEYb2yu5j8fV7BsVz1aw7TiDK6cXMglY/NOWFdTJBTj4zf3s/bNMpSCSRcOYuL5RVhs5kNl9q3/mNcf+iNBn5ezrv88Ey+e36ZLp7I5wPWPrqDBG+aRmyYzc2hWm/fQMYP6BVsI7Woi65bROHI88ZZD4x6Y/yeYcH2X6vpO2Tv8aOmPsJvt/OrMXzG74NihIoQ4tUgwdNGB5gD/XVvJwjUV7K33kWS38NmJ+Vw3tYjR+cc+y6crPPUBPnx+N7s/riUp3c7UecUMn553qBvL72nhjb8+wJ41HzF4wmQu+tq3cad90jqobgly8+Mr2Vfv575rxvOZcW2vczBCUer+uoFoQ5DsL4/FlhWD526GvUvgzO/B2T+Gdk6RPdKelj3c8d4d7Grexa1jbuX2ibdjMZ348RghRN+QYOgmrTVr9jfxr4/KeGVDFaGowfjCVK46YyCjhmVQGY1RE47QEI7SHI0R05pY4nNMsphJtZhJsZgZYLMy0GljoMNGusXc5q//Azub+OA/u6nd5yEj3830z5YweGwmSim01qx/6zWWPPEYFrudcz7/FUbMnnvo9S3+CF96YhWr9zdxz7zR3DJzcJv6xzwhav+6ASMQJfsr47Dl2ODl78DaJ2HUZXDZX455IRxAMBrkNx/9hv/s/A+Tcibx27N+K11LQpwmJBiOkz9m8OqBRp5ZU876LfVEPGG0WRHLcxIb6Eal2Ei1mDErhTnxO781ZuCPHb0EZ7LZxAi3k1FJDkYnORmT7GS020HF+gZWvLiH5ho/eUNTmXH5UPJK4q2Thspy3vjrA1Tt2MaQyVM5/0u3kZQRvxAtGInxjafX8taWGm47u4Q7zh/eZvA82hik7q/r0YYm+yvjsGY7YfmD8NZPIXNYfNwha1iXPodX9rzCvcvvxWa2cc+Mezh30NFnTwkhTi0SDN2gtea9xlaeq27kjQYP/piB06SYkuIiPxj/S3/DjgZCUYMxBSncMG0Q88bnk3TYWETE0HiiMapCYcqD8dueQJit3gBbvAFaE8HhNCkmprg5I8XFgIoQodcriTWFKR6fxfTPlpCR58YwYqx97SWWPfMkZquFuTd/mdFzzkUpRTRm8JMXN/P0R2VcMnYAv79qPC7bYfWo81P3tw0okyL7f8ZhyXTCniWw8AsQDcPlD8PIeV36XPa17OP773+frY1bmV8yn7um3kWy7ditDiHEyemkCQal1EXAA4AZeExr/Zsjnv8u8CUgCtQBX9Rad7zuJScuGCKGZmFNI38tr2O7L0iG1cxnstOYn5PGGalu7If1y7cEIry4rpJ/rSxjW3UrbpuZ+RMKuH5qUbvXGRxOa015MMz61gCrWnysbPGyyRsgpsGqFCOiZnK3eymuCDFnVDZTLy4mJctJU1Ulb/z1T1Ru20zRmHGc88WvkVkwEK01jy3dy69f28rIvBQevXkK+WmfzH0UrvJR/+gGlMVE1pfHYs12QUtFfNyhcg3M+jacczeYrcf+jGIR/rrhrzy28TFyXbn8cvYvOWPAGcf9mQsh+s5JEQxKKTOwAzgfqABWAddprbccVuZsYKXW2q+U+howV2t9TWf7/bTBoLXmrQYPv9h9gJ3+EKPcDr5alMNnc9KwHWOQVmvNuvJm/rWyjJc2HCAYMRhbkMp1U4uYP6FtK6IzvmiMNR4/7ze1srjRw2ZvfGEfd9BgWHWE81KSuGVOMVkZDta/9RrLnnmCSCjElHmXM/2Ka7DaHSzeVss3n16L3WrmbzdNZvKgTwasw1U+6v9vIwBZt47FlueGaAhe+wGs+Xt8Ar7PPQoZQ7pU3/V16/nxsh9T5injplE3cfvE22UiPiFOMSdLMMwA7tFaX5h4/EMArfWvOyg/EXhQaz2rs/1+mmCoDUX4wY4KXqtvYajLzt1D8rkwK+W4rvr1BCO8sPaTVoTLZuayCflcPamAkgwrgUCAUCjUZh0Ei8WCw+E4dDObzYfq9V5TK29VNfNuowefCaxRzcSomatLB3Buqpl1zy5gy/vvkpKdw9mf/x+GTpnGzppWvvTEag40B7j70lHcPGPQoWOJ1Pmpf2wjRtgg+4tjsA1MdANt/i+89C0wYnDpH2D8tV06Xn/Ez31r7uPZ7c9SmFTIT2f8VKbTEOIUcrIEw5XARVrrLyUe3wRM01rf3kH5B4FqrfX/tvPcV4CvABQVFU3ev7/T3qZ2vVbXzHe3leM3DL5fnMdXCrOxdnTls78R6ndA3TZo2g/eGmitBm8thFshEkBH/PF7rYhiIqJNBLWVVu3Ei5MYVlAmAspFC8mHbh6S8eJGo3C73aSmph66ZWVlkZ6VxeaInX9vr2e5KYLXacKsYUaKm7kqgOnffyewdyeDJ0xmzo1fxJKZxx3PreedbbVcOjaP33xuLMmOeDdRtDFI3WMbMbxhMq4fiXNEYmK+5nJ4/itQ9iGMvQou+X/gTG//szjCqupV/Hz5z9nv2c/8kvncOeVO0hxp3f4+hBC965QLBqXUjcDtwBytdaiz/R5vi+HN+hYe2F/D/SOKGOY+bOqHWBSq1kHZCihfAeWrwFt96GmtzGh7OtqaTtScTEvERHNE0xTWhDGj0NjNJlwOGza7lUA4QsDnxRbzkq68DDC14NCBNnWJme34nYW0WHOpVxlURZLZH3BRE01CJ9ZSSkpKIj09i3Ijg7U6mR15yTQlWbAqmBDxk7/8bQbv3MCUM+cy/crreWpjE//vje0MTHfylxsmHboOI+YJU79gM5EDXtIuG0rS9Lx4JYwYLL0P3vs1uLPhM3+EEZd06bMMxUL8bf3f+Pumv5NiT+HOM+7k0uJLZb4lIU5iJ0swdKkrSSl1HvBn4qHQ8VqYCZ+mK8nQGpNSEPLCjtdhxxuw6y0INAEQs2QRCmcSqNUEyrwEa6KE/WaqBuSzp2QIB/Lz0SYTLp+P/AMHyKmpJbOxAac/gIofDOaMDKyFhewcOIqXU4bxZjAJmxHgrOwAnys1MSPTj8OzD+q2x1slLeWH6qetLoJpw2hyDKKKHHb5k9jeAAYKDTS4ctiTNYjdhbm0WC1YDYMh+7Yyev9WrpkwDvuYM/nWws00+SJ878JSbp09BLNJYYRiND69jeC2RpLOKiT1osGfTPBXtR5euA1qNsKYK+Hi34G7a2s17Gjawc8//Dkb6jcwKWcSP5z2Q0ZkjDiu70YI0bNOlmCwEB98PheoJD74fL3WevNhZSYCC4m3LHZ2Zb/HHQyGAfuWwvqnYcsiiPgwzMkEvNm0bAngKzeIBs2Y09OxjxiOtXgIe1NT+DgQoCEYxGWzM3ZIMWOGDSM7JQUiEXQ4jBEIEGtqItrYRKyxkWhdLeHyCiJlZUSqq/GZbCweOIlXB09nb2o+Th3louQg14zPZcqZE7C4rFC/E2q3xH9JH1gH1Rsg4gfiYRHKHE2ds4RdwWzW1lppiWmqUzLYm1vEnpwCvBYr1kiIkQf2MC83iw21yby1tYGpxRn84arxDMxwoWOa5pd241tRhWNUJhlXl2JyJAbLo2H44H5Y8rv4ug4X/xbGfA660AKIGTFe2PUCD3z8AC3hFq4qvYrbJ9wu3UtCnGROimBIVOQS4H7ip6s+rrX+pVLqXmC11nqRUuptYCxQlXhJmdZ6fmf7PO5gePeX8P7v0BYXPm8BDSs8+GtMWPMKcM+ejXv6NJyTp2DJyWbv3r289tpr1NXVkZ2dzZlnnsno0aMPDRZ3lQ6HCVdUEtq1k+DW7Xy8q5oX/CksTi8lZLFR0lzBPN9uLh3oIGviWFwTJmArKYm3Dxp2xUOicg2Ur4TqjaBjaBTR9FLKwkVsaclgpymTLRkF7MopZG9WPiGrDWc4yNhggO27DCyeKD+bN4qrJhcC4P3gAC2v7sGS6STzplFYcw6bQLBmM7x4GxxYC8Vz4JLfQ3Zpl47VE/bw8LqHeXrb07itbm6bcBtXDb8Kq+nYp8UKIXreSRMMPeF4g8H7yj/xPfsgTWuaMaVkkDp/HimXXIJj3LhDfeNer5dXX32VLVu2kJaWxvnnn8/IkSMxdWGuoe5ormti4VvreHZrMzsjNhzRMHMq1nJ+2SrGGs24p07FNX0a7mnTsA0ZEq9fyAuVq6Fs5SdjIeHW+P6MPPZFhrLfUcS76UNZk1vKvswBRM0WXMEgpio/oy1uHrhkNIOzkgjtaabhX9vQYYP0q0pxjT1scj4jFj+l9Z17IeyHGbfBWXd2aflQgF1Nu/jNR79hZfVKBqUM4hsTv8EFgy6Q8Qch+pgEQztaXnmFhkcfI+Omm0j5zKWY7G3XQNi+fTsvvvgioVCIs846i5kzZ2K19uxfu1pr1le08K+V+3lpXSWBqGYgAc4tX8PZW94jJ9CMOTsL99Rp8aCYMQNbYfwvf4wY1GyCfcswdr+H3rsMcyze/VRvDKTCPojFjhJeKpjDtuxiDJMJl8/PRKW4Z8YoRpodND61lXBZK+4ZeaRdUoyyHtYi8tXD2z+Dtf+ElAI496cw9uouTcintWZJxRIe+PgBdjXvYkzmGL4z+TtMzZvaEx+jEKILJBjaoQ0DlDrqL9dYLMZbb73FihUryM3N5YorriA3N/dEVbfLvKEor22sYuGaClbubUQBU5NjXNi0nSmrXsVaEz9TyjZoEO4zz8Q9exbuqVMxuRJdQbEIunItzStfJ7p9MemRTVhUGANFlc5huWM0Lxeew1v504mZrAwI+rk4ycYV9TbyVrZgyXWRed0IrAPcbStW/hG8emf8zK0B4+D8e6Hk7C4dU8yI8dKel/jLur9Q7atmVv4svjHpG4zOHH3iPjghRJdIMHRRIBBg4cKF7N69m6lTp3LBBRdgsfT9VNPljX7+83EF//m4gvLGAEl2CxcNdnNxqJyha98jsPIjdDCIslpxTplM0uwzcc+ejb102KHga9hfT9kbr2PsWkKBZS25tl0oNAHDynr7cF7PPZNXCuZQ7sijyNvCmVVe5lWZmHTuWJJnFrRdltQwYNN/4t1LLWUw9Dw47+cwYEyXjicUC/HMtmd4ZMMjeMIe5hTO4X/G/Q9js8f2xMcnhGiHBEMXeDwennjiCRobG7n00kuZPHnyCajdiWUYmo/2NfKfNRW8srEKfzjGoEwXl43N5QJTA9nrVuBbupTQzvgJXZacHNxnziZp9mzcM2diTk0l0Bpm89JKdi7dTrp/FYNdKyi0riXJ5AOg3JzDksxpvJEzk+Up48ny+Jjd7OeasYMZP3EUNpvtkwpFQ/DRo/D+/4NgC4y5As76PuR07RTV1nArz2x7hgVbFtASamFW/iy+Ov6rTMiZcKI/OiHEESQYjqG5uZkFCxbg8/m47rrrKC4uPuZrtDaIRBqJRFqIRJqIRFswjBDaiKJ1FK1jKGXBZLZjMtkxmxyYzU6s1nSs1gwsluObeuMgXyjK65viK88t39OA1jCmIIX54/O5KM9KyoZVeJcuw/fhhxgeD5hMOMePjwfFmWdhGzGCsq3NbH6/kv2b60lTmylOfpcBxgYGOpuwmgwiyszHyWN4K3May1In4g2lMNtiMH9gHqOHl5KampgwMNAEHzwQD4mwD0ZfDnO+Dzkju3YsER/Pbn+WBZsX0BhsZNqAaXx+zOeZlT9LBqmF6CESDJ1oamri73//O+FwmBtvvJHCg4O5CYYRwefbiad1I17vVgKBssStEq3Dx/2+SlmwWtOx2bJxOPITtwIcjnycjkKczsFYrSld2leNJ8jLG6pYtK6S9RUtKAVnDM7gsgn5XDwyG8fu7fiWLsO7dCnBTZtAa8zp6bhnzybprDOJjZzC9k1+tn5wAL/Hi1ltJD32DiWOCoqSWsi1x892arCk8V7GFD5IncS+WAEl2LmgMJeRpaXk5+djCjTF13z46JFEQHwWZn8X8sZ16Tj8ET8LdyxkweYF1AZqGZo2lJtH3cylQy7FZrYdewdCiC6TYOiA1+vl8ccfx+/3c8stt5CXl4fWMTytm2hsWEpj0wd4POswjHgAmM1JuFyDcDqKcDoHYnfkYbWkJVoBqZhMDkwmK/GJZM1oHcUwghhGiJgRIhbzEQk3EYk0Eo40EQk3EArXEQxWEgweIBbztqmfzZaFy1WC2zUElzt+73YPw27P6/Av6X31PhatP8CL6yrZXefDYlKcVZrNvPF5nDsyF5e/Fd8HH+BduhTfsg+INTaCUjjGjME5+0zq86ey54CNsq2NGJG9BCJrSYntZJC7iTG5DnIs+3HoFgC2uov5IGUiuxlEJJbNzPwCRg4fzpC8dOxrHoGVf4OwF4bMhZnfgJJzu3SRXCQW4fV9r7Ng8wK2N20n05HJtSOu5Zrh15Du6NocTkKIzkkwtCMUCrFgwQJqa2u56aabSEtrpLrmJWprXyUcrgMgOXk06WnTSU4ZS0ryWJzOIpQ6sdcwHKS1Jhr1JEKiAr9/Hz7/Hvz+3fh8u4lGWw6VtVhSSUoaQXLyqPh90ijc7qGYTLY2+9tS5WHRugO8tP4AB1qCWM2KWUOzuHjMAM4fNYB0p4Xg5i34li3F+/5SAuvXg2FgSk3FNP1sagaeyd6mFJpr64iENxEOb8KpQ0zNHc7QlAAR+3qyYluw6TAhZWVVyhg2WYdRE82jKGMoo4YUURpYQ/r6R+LzTeWMghm3w9grwWJv72M46jNZWb2SBZsXsKxyGTaTjQsHX8jVw69mfPZ46WYS4lOQYGjHO++8w/Lli7n44jTCkTfw+/dgMtnIzJxLTvZFZGTMxmbr2hxBPU1rTSTSiM+/B593B17vVlq9W/F6t2EY8bUblLLgdg8lOWkUySnjSEkZR3LSCEwmO4ahWVfRzOubqnltUxXljQFMCqYVZ3Lx2AFcOHoAuSkOYi0t+D78MD42sXQp0bo6NBAYdzbVg+eyz5tOLLSfUGgjybqRCRlzGeguIurYiL9wJ8meVRR64gPf9dY0PnaOYq8xEGUfwpAMJ6V1b1DY/CEmdzZMuhkm3wJpRV36DHY37+bpbU/z8p6X8UV8lKaXcnXp1Xym5DO4re5j70AI0YYEQzvKyp9m165fobWflJTxFORfT07OhVgsp85ylVrH8Pv34/VuodW7Da93Cx7PJiKRBgCUspKUNIKUlHGkJI8jJWUsLlcJW6p8h0Jid50PpWBSUToXjR7AuSNzGJKdhNaa0Pbt8S6n95fiX7uWmKFoKJhK5eA5NKgUjPAO0lUz41PHkenIJ2DyERvpIzSwnsje9xl84AOyIvEJCXfbC9lmLqEpNoAkq51i7ypKKCOp9EyYcisMPRdMx55ixB/x88reV3hu+3Nsa9yGy+LikiGXcFnJZdKKEKIbJBja0dj4AVVV/6Gw8GZSUyec+Ir1Ea01oVAVLZ71tHo24PFswNO66dD4hdnsIjlpdDwsUsZRGxrOuzs0r2+uYWuVB4DiLDfnjsjhnJE5nDE4A6vZRMzrxbd8+aFB7EBdC3VZ46konEWLI4lcUz2jkwpIs2XSEqmjKbORtNlDCWVaqN21mLT9S5jYtB6HDhNWFrbaSthvFBCNOkjVTQxxhxg45WIsE6/tUitCa82m+k08u/1Z3tj3BsFYkEEpg5g3ZB7zSuaRn5Tfo5+zEKc6CYZ+TmsDv39vIiQ24PFsxOvdfGhQ3WrNICVlPH41mfU1Q/lwv4uVez2EYwbJDgtzSrM5d2QOc0tzSHfb0FoT3r0b7/tL8S1bSuO67dSlj+VA9mRSMjMZ7nCSYk3CG2lme+tqKLExdPoMAsNK2btnBWrPYkbWLGe0bzcAfuVgp2UQtdEsYtpMRnomJZPOInPqVWA/dgvOF/Hx5r43WbR7Eatr4v8Wpg6YyrySeZxXdB5Jtq7N6yREfyLBII5iGBG8vu3xsPCsx+NZj8+3C4j/e1DWoez2zmFtbSkry9w0+AxMiS6ns0qzmT0si3EFqVjMJgy/H9/KlbS89z4VH66jNZZNrGg2+Zl5ZFgdBGN+dnrWsNe7kZxRwyidNpOksRNZ52uhesd7JJd/yLSmjyn1x1fi8ykHe9VAGox0VEoBmePOZ/DUS3GnpB3zuCpaK3h5z8ss2r2I8tZybCYbMwtmcuHgC5lbOFdCQogECQbRJdFoK57WTW3CIhSqxtCKstbBbGk+i/V1I9jd4EKjSHFYmDU0izOHZXPmsCwGZsTnaQpXV7Pz9ffY8/Zy3L5UsgrPIMOdhqENKvy72dmykvpQJZmFgxl6xlQKJkxhf3YBH1XupXX3EobVrmJW81pKAhUAeJWTfRTSbB6AKphE5hmfY1BJKU6ns8Nj0Vqzvm49b+x7gzf3v0mtvxabycasgllcMPgCCQnR70kwiOMWDFXT6tlAi2cDHs86PJ6NNAc0WxtL2dI4li0NI2kIxM8KGpRh56zSAcwelsW04gzSXDaCkSjvv7eO/S8vp9SfwqC0IixmC80RL3u8W9jvWUHYCGBzuCiedAYlU6aiSkfzURg2VO5C7VvKxMaPmdW8lkHB+DIdAexUMoBm5yAoPou0sRcycMhw7Pb2T4E1tMGGug1tQsJqsjJ1wFTmDJzD3MK55CXl9dpnKsTJQIJBnDDx8Yo9eDzrafGsp6VlPbtqGtlUP5TNDSPY3lRKKGZDoRmWY2FmyQCmlwxgWnEGkZjBax9XEHp3K1NazeTYkzG0wYGQh73ebVR5P0In1sNOy8qjZOpUCsdNorFwCB/4Qmwu34H9wEdMadnItOYNjPDvxYQmhokqcmhylRAtnI575DnkD5tAUtLRLQJDG6yvW8/b+99mScUS9nvi3Vel6aXMKZzD3IFzGZM1BlMPXa8ixMlCgkH0qFgslDhVdj31TRv4eH8dG6uT2dY0lN3NxUSM+IV3JZkxphWnMat0CIOzktn0cRWhdXVM8cRIMZkJxiIcCHooC1RQF9iEET0AGChlIju3gCHTppMz8Qz2ZObxQVMrGw7swd2wjTM88aCY2LoVh44AUE8a9bZBhHPGYSueSdaYuWTmHH3F+N6WvSwpX8J7Fe+xrnYdMR0jw5HB9LzpzMifwYy8GeS6e3/adSF6mgSD6HWRSBMez0bqmtbz8b4KPi6LsKW+gF3NxYSNeJdPYWqQ8QU2SnNycNTZGLjXywi/xobCG4tQH/JSEWikOrQPI1qGjtUCYFImMjNyGDRmLBnTZ1MxYCCrGmpZXVeD9hxgimcz01o2cEbLJjITV4xHMFNrGoA3dTgUTMZdOofs0jOwOxyH6twSamFp5VKWVixlRdUKGoONAAxJHXIoJKYMmCIX1InTggSD6HNaa4LBShqa1rF6zy4+2tfKlhoHu5qL8EbiXT4ua5jSzACDbQZDWm3MacokExtNOkZtLIAv5OWAr5ZgtBojWomO1XHwLKoUZyoDhwwhZ/pM6oaUsLbhAMs9QWrDUcb4djLRs42pLRsZ7d2FU4cA8OGk3jaQQMZozAMnkzL8TLIGj8FssWBog51NO1lRtYLlB5azpmYNwVgQi7IwJmsMk3InMTl3MhNyJpBi69pkh0KcTCQYxEnJMKJ4vTvYUr6R1ftqWF8ZYVtdChWtA9CYUBjkuZsYbPcywtCMiSUx1JdDmeHAZ4piBD14PLUEIvXoaCVGtAqIAWA32chKzSCntBT/qBFUuGCdP8BanYI76mOSZwuTPVs4o2UzQ4LlmBIB48VFs2Mg4cyRmAsnkTxsFmnFE4gQY13tOpYfWM7qmtVsrt9MVEdRKErTS5mcO/lQWGQ5szo5aiFODhIM4pShdYy6pj18tHs7q/fVsakqyq6GZJpD8QvdFAZ5rnpK7B5KTVEKDTuEUwmRjz0UItZYhc9fjxGtRceq0UZi8kENLrODnLR0nEMGUl+UzR63jXXazQ5bLiWBCsZ5dzDes53xrTsoCZRhSYRMEBtN9kKCaaWQNwH3kDNwDhnP1ta9rKldw5qaNWyo20AgGh84L0wqZGzWWMZkjWFc9jhGZIzAYXG0e7xC9BUJBnHKK6urZM2ebawrq2FLdYid9W6aQ/G+foXBAFcdRc5G8i0hkoFkUzLJ5GOvsRBprCEQqMOI1mDEqiFx5hMaXMpOmttFtDCXxvx0KlNdbHNmsd1VRFGwhrHeHYxv3c4kzzaG+ffhSKzBYaBoMWfhSxqEkTUCS8FY6jNzWWtqYn3DJjbVb6LKFz+91qIsDEsfxpisMYcCozi1GIup75eNFf2XBIM4LR1o8rBmzw7WlVWxtaKV3Y1magKfnKLqMAcpTKoi395CuiVMus1KjimZlAYrvvIwAZ8HHa3DiNWB/mQtDKthJtlsxchNpzE/m6rMVHZnDGBrUjEDwo2M8O1hhG8vY727GOHdQ0G45lBXVAQLLdZcAilDCKYXU5eawXanwYrYATY2bKU1El/0yG62MzRtKCMyRlCaXnroXi66E71FgkH0G95QlK0VtazbupNN5TXsaIyw3+/EH/vk4rdUWwt57hqyHS1k2EJkWi3kRsy46qN4qxVBTwAdrUMbjRwc3EaD3bAQS0nDl51KU1Ya1VmZ7MsYSJUzmyHBSkb69jDSt4exrbso9e0lO9p06D2jmGm1ZuN1DaDBlUml08kma4glsQOUx1oPlStMKoyHREYpQ9OGUpJawsCUgVhN1t76CEU/IcEg+jWtNRWVrWzYVM7GfWXsbGqhPGhQEXHiNz5Z3MhuDjHAXcMAVx1Zdj9ZVoOMcJTk1iBGkybUEiXUEkDHWjg8MEzYCKRm0pqVRmNWOtXZA6hMz8Nnc1ISqGBooIyh/jJGePdR4i9jYKj60PgFgE+5abFlU29Pp8JmY7s5xBrVyB6rQavFhEVZKEopoiSthOLUYkpSSyhJK2FQyiAZuxDHTYJBiCNorYk2BDiwu4kdexrZXFnFXl8zZUaU8piN+mjbeZiSrF5yXHVkOxvItPvJUCHSI0GS/H4snhBhT4xQc4hY2ADiF9HFTFa8SVl40zLxpiXTlJFObUYOVak5ZMVaKPGXMSxQxjBfGaW+/RQHysmIedq8r185abakUmdxU242s90UYY/NYJcjRo3FRF5yIUUpRRQlFx26H5gykMKkQlknW3RKgkGILtBaE/OEidb4aalsZef+ZnZVt1Dpa6ZCBTiAwQHDQn3UjuaTKTPs5hDZznqynI1k2FtINQVI1UGSIwFcIR92r49wa4SI10TEZ8UwFD5nMi2pObSmZuJJSaM1NYWGtAx8bjc5sSYGBw4wKHiAwYFKhvgrGBQ4QH64rk1LI4qZFlMSDSYn1SYLlWZFhVVRbtOU2QyMtHxy0gbHQyO5iPykfPKT8slz55HhyJBFjfo5CQYhPgWtNbGWMNEaH5EaP74aH7sPeChv9FMV9lNu8VOpotRoRW3Mht9oOx5gUVEyHE3xm7ORdLOPVOUjKRrAGQ7gDLZi8fmJBsyEfVY8RipNrkw8yVl4kzPwpKTSkpJGU2oabkuYokRgDApWURCqoShQTUGwhpxIA2ba/n/2KQeNykmtyUaN2Uq1xUyVRVHrMBNOzcCWXkR6ykDykwvIc+fFw8OdT7YrW86aOs1JMAjRQ3QkRrQhSLQhQLQ+SLjeT12VlwMNPuoCEWpUhAPmANWmMDVaU2NYaIrZ2rQ4ACymCGn2lsTNQ6qllWSTjyTtxx314Yq04gz6iAXNNJJFo8qkyZJBqyMDb1IarUmp+JKSsNuj5MQaKQjWUBCqpTBUQ0GwhoHBavJDtaTE/EcdQxArzSYnzVhpMFupN1uptZjwOFwEk9KIpediyygkNbmAXPcAsp3Z5LhyyHHlkGZPk5bHKUqCQYg+oGMGsZYwseYQ0eYgsaYQocYArXV+apoCNHrDNOgYdaYQteYg9SpKPZo6baIpZiVgHH0mks0UJtnWSorNS7KtlWSrD7fy41J+XNqPM+ZHGZogLgImFy3WDJpNabRa02m1pYLDgssSIlX5yQ03MCBcT26oIf5zqJ4BoTpyww04E9drHC6MGY/JhUfZ8CgLzcpEs8WK3+Yk6EwhmpQGyVlY0vNwpReRmjSADEfGoVuaPQ1zF9b1Fr2js2CQtqIQPUSZTVgyHFgyHNhJPbQ9GxhCvIvK8EXi4dEaxvCEiXpCBJuCBJqCeFuCNHrDNIeiNGBQawrRYIZGw05TyEZTMJMKw4wnZiOq2/+F67L4SbZ5SbZ5ybbW4zIFcJmC2AiB0lQqC/tN2YTMhQRS3PhtbrzWVJTVSpIOkYKfzFgLueEGMsPNZEWayYo0kRlpZmi4iaxgA45AGFqOfu+AsuFRLlpNNsqUhS3KRMBsIWh1ELa7iTlT0EnpmJIysaXm4EjNw51SQHIiRFLtqSTbkqVLqw/IJy5EH1FKYU6yYU5qe/ZQ6hHltJEIEE8iQLwRDF+EkCdEoCVEqCVIqzdMSyCCNxyjxYjRpCI0mkM0mqAp5qQp4KTVn02lYcZrWAm20xo5yKRiuCwBQtYIQWuEVquJClMGypQOJhOGyUTUZSWSbCNisqEsZlymCC5ThGQVJBU/WdHmQ0GSEWkhLdpKQdRDWqSRtEDroQsC2+NXdlpNThqVnTKTDb+yEjJZCFusRC02DJsTw5GEcqaiXClY3BnYk7NxpObgdOfgcmWRZEsmyZaE2+KWVspxkGAQ4iSnTApzsg1zctsASe6gvBGOYfji4RHzReJBEogS9oYJtIYJecMEfGFafRFaQ1H84Rj+aIxWFabZHKLFFMGjNJ6YA0/MQWvAhFeb8Btm/IaVWAetk8PZzCHsljA2SwSLJYbZEkNZFNqm0CaFU0VIUiFcphBJKoybEG4VJE37SDO8pBmtpBmtpEc9ZBh+kqJNpIR9JLczRnKkGCZ8Jgd+ZafOZMdvshFUVsKm+C1qthIzW9FWO1jtKJsDZXOjHEmYnalYnCnY3RnYkzJwuDNxurJw2lNwWV04Lc5+0YI5/Y9QiH7GZDNjspkhvesXv2mt0REDHYhiBKIY/sR9MEokECXoDRPyRwj5wvj9UVoDYXyhKP6wgT8SJahDBAgTUBG8KoJXGXgx8MUUvpgFX9BGQJsIaRNBw0y9YSGiu/7rR6GxWKJYzVGSTAFSzH6SVJAkFcBlCuNWQZII4lIhkgiSjJ8kHSAZP8mGnyQjQHrMi0sHcRsB3EYAC0aX3z+sLASVjUZlJ2iyEcZCWMVvEWUhYjITNVmJmCxEzRYMsxXDbENbbCizDWWxoax2rDYHZpsLsz0JiyN+szpSsDlTcbjScDpTsdtTcNjc2My2PrvivVeDQSl1EfAAYAYe01r/5ojn7cATwGSgAbhGa72vN+soRH+klELZzGAzY05tf+3srtAxjQ5FMcIxdDCGEYqhw/H7cDBK0B8hHIgSCkQIBqJ4AxH8oSiBUIxgJEYgEiZohAnpECEjQogIQWKEMAgqgwCaIJqA4SCsndSTTjgROGFtImyYCWszYcNyjJaNxkYUF0FchHCpIO5EsLhI/GxKbFchXKYgTh2K/2yEcKgwDsLYieAkTDphHDq+zUEIhw53K3iOFMFMRFnwKWs8eNrczESxEFVmdqadwdVffei436cjvRYMSikz8BfgfKACWKWUWqS13nJYsVuBJq31UKXUtcBvgWt6q45CiE9HmRXKZcXkOvovXden3LeOGvFbxDjq50goRigYJRyKEg5GiYYNQsFoInQiBMMxguEYgXCEcCxMMBYiFIsQiSnChoWw4SBMjIiOESNGFAMfmhatiQJRIKIhiiKKIoyKBxEmIrrtLaZNRA0zShtYjQh2YjhUGDsHgyN8KFgcRNo8thHFRgSbimAjip1I4nF8uz1xO/h4dw91a/Vmi2EqsEtrvQdAKfUMcBlweDBcBtyT+Hkh8KBSSulT/ZxaIcSnpiwmlMUE7fSQ2YHempdWaw2GRsc0RI14KymmIWYQixhEozEiYYNYJEYkEiMcihEIxQiFowRDMULRGKFIhHAkTCgaJhKLEolGErcokViMaCxKqxGl2YgSM2LEjCgxHcMwDGLaIIZBTGvGZ2X3yDH2ZjAUAOWHPa4ApnVURmsdVUq1AJlA/eGFlFJfAb4CUFRU1FP1FUKIoyilwKxQZsDWtrvKQjykTnWmYxc5+WitH9FaT9FaT8nO7pnEFEKI/qo3g6ESGHjY48LEtnbLKKUsxE/pbuiV2gkhhAB6NxhWAcOUUsVKKRtwLbDoiDKLgFsSP18JvCvjC0II0bt6bYwhMWZwO/AG8dNVH9dab1ZK3Qus1lovAv4PeFIptQtoJB4eQgghelGvXsegtX4VePWIbT897OcgcFVv1kkIIURbp+TgsxBCiJ4jwSCEEKINCQYhhBBtnPIL9Sil6oD9x/nyLI64eK4fkGPuH+SY+4dPc8yDtNbtXgh2ygfDp6GUWt3RCkanKznm/kGOuX/oqWOWriQhhBBtSDAIIYRoo78HwyN9XYE+IMfcP8gx9w89csz9eoxBCCHE0fp7i0EIIcQRJBiEEEK00S+DQSl1kVJqu1Jql1Lqrr6uz4milBqolFqslNqilNqslPpWYnuGUuotpdTOxH16YrtSSv0p8TlsUEpN6tsjOH5KKbNSaq1S6uXE42Kl1MrEsT2bmNEXpZQ98XhX4vnBfVrx46SUSlNKLVRKbVNKbVVKzTjdv2el1HcS/643KaWeVko5TrfvWSn1uFKqVim16bBt3f5elVK3JMrvVErd0t57dabfBcNha09fDIwCrlNKjerbWp0wUeAOrfUoYDpwW+LY7gLe0VoPA95JPIb4ZzAscfsK8HDvV/mE+Raw9bDHvwX+qLUeCjQRX08cDltXHPhjotyp6AHgda31CGA88WM/bb9npVQB8E1gitZ6DPEZmg+uC386fc//AC46Ylu3vlelVAbwM+IrZE4FfnYwTLpMa92vbsAM4I3DHv8Q+GFf16uHjvVF4HxgO5CX2JYHbE/8/DfgusPKHyp3Kt2IL/r0DnAO8DKgiF8NajnyOyc+7fuMxM+WRDnV18fQzeNNBfYeWe/T+Xvmk2V/MxLf28vAhafj9wwMBjYd7/cKXAf87bDtbcp15dbvWgy0v/Z0QR/Vpcckms4TgZVArta6KvFUNZCb+Pl0+SzuB74PGInHmUCz1jqaeHz4cbVZVxw4uK74qaQYqAP+nug+e0wp5eY0/p611pXA74EyoIr497aG0/t7Pqi73+un/r77YzCc9pRSScB/gG9rrT2HP6fjf0KcNucoK6U+A9Rqrdf0dV16kQWYBDystZ4I+PikewE4Lb/ndOAy4qGYD7g5usvltNdb32t/DIaurD19ylJKWYmHwlNa6+cTm2uUUnmJ5/OA2sT20+GzmAXMV0rtA54h3p30AJCWWDcc2h7X6bCueAVQobVemXi8kHhQnM7f83nAXq11ndY6AjxP/Ls/nb/ng7r7vX7q77s/BkNX1p4+JSmlFPHlUbdqre877KnD19K+hfjYw8HtNyfObpgOtBzWZD0laK1/qLUu1FoPJv5dvqu1vgFYTHzdcDj6mE/pdcW11tVAuVJqeGLTucAWTuPvmXgX0nSllCvx7/zgMZ+23/Nhuvu9vgFcoJRKT7S0Lkhs67q+Hmjpo8GdS4AdwG7gx31dnxN4XLOJNzM3AOsSt0uI962+A+wE3gYyEuUV8TO0dgMbiZ/x0efH8SmOfy7wcuLnIcBHwC7g34A9sd2ReLwr8fyQvq73cR7rBGB14rt+AUg/3b9n4OfANmAT8CRgP92+Z+Bp4mMoEeItw1uP53sFvpg49l3AF7pbD5kSQwghRBv9sStJCCFEJyQYhBBCtCHBIIQQog0JBiGEEG1IMAghhGhDgkEIIUQbEgyiX1JKZSql1iVu1UqpysTPXqXUQz3wfv9QSu1VSn31BOzr/yXq/L0TUTchjmQ5dhEhTj9a6wbiF4mhlLoH8Gqtf9/Db3un1nrhp92J1vpOpZTvRFRIiPZIi0GIwyil5qpPFvu5Rym1QCm1VCm1Xyl1hVLqd0qpjUqp1xPzUqGUmqyUWqKUWqOUeuPgvDbHeJ9/KKUeVkqtUErtSbzv4yq+6M4/EmXMiXKbEu/5nR49eCESJBiE6FwJ8Yn55gP/BBZrrccCAeDSRDj8GbhSaz0ZeBz4ZRf3nU58DYHvEJ/35o/AaGCsUmoC8RZNgdZ6TOI9/36iDkqIzkhXkhCde01rHVFKbSS+atjrie0biS+oMhwYA7wVn9sNM/G5brriJa21Tuy7Rmu9EUAptTmx7yXAEKXUn4FXgDdPyBEJcQwSDEJ0LgSgtTaUUhH9yeRiBvH/PwrYrLWecbz7TuwrdNh2g/iqZE1KqfHEVyr7KnA18cnRhOhR0pUkxKezHchWSs2A+HoYSqnRJ2LHSqkswKS1/g9wN/E1F4TocdJiEOJT0FqHlVJXAn9SSqUS/z91P7D5BOy+gPjynQf/gPvhCdinEMck024L0QsSZxq9fCJOV03s7x565xRb0Q9JV5IQvaMF+MWJusANuJH4Ws9CnHDSYhBCCNGGtBiEEEK0IcEghBCiDQkGIYQQbUgwCCGEaOP/A104knOTaA02AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(model.t, model.output.T);\n",
    "plt.xlabel(\"Time [ms]\")\n",
    "plt.ylabel(\"Activity $x$\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## BOLD and autochunk"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we've followed the model implementation guidelines, the model is also compatible with chunkwise integration and can produce a BOLD signal. Let's try it out!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {},
   "outputs": [],
   "source": [
    "model.params.duration=200000\n",
    "model.run(chunkwise=True, append_outputs=True, bold=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'BOLD activity')"
      ]
     },
     "execution_count": 75,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZ0AAAEGCAYAAAC+fkgiAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/Il7ecAAAACXBIWXMAAAsTAAALEwEAmpwYAAAnsklEQVR4nO3de7xcZX3v8c93z2XnQq4QIBAwQSK+ACuXHMTj5bRaETwqiqjh2IqVilTxWKz1hceKFmtb7LG2Cl7oEVGPCoha04riBcXqESSxyE2RcFESA4SEJIQk+zLzO3+sZ3Ymm9kzk71n1sw23/frNa+95lnPWus3a+/ML8+znrUeRQRmZmZ5GOh1AGZmtu9w0jEzs9w46ZiZWW6cdMzMLDdOOmZmlptirwPoBwcccEAsXbq012GYmU0ra9aseTQiFu3NNk46wNKlS1m9enWvwzAzm1Yk/Xpvt3H3mpmZ5cZJx8zMcuOkY2ZmuXHSMTOz3DjpmJlZbpx0zMwsN046ZmaWGyedDtnwq3v41qeu7HUYZmZ9zUmnQ955zTe4+AHx8X/8SK9DMTPrW046HfKc4q18v/wOvrdN/OaX9/Q6HDOzvuSk0yFzYgeDGuVD1U9ywdXfIqrVXodkZtZ3nHQ6pECFaohlAw/x6pHv8t6//rteh2Rm1necdDqkEFUeZxY/mf1cVhZ/wLbYwaMbH+11WGZmfcVJp0MKUWGEAie+5UvcVXwKH9AV3PTTH/Y6LDOzvuKk0yEDVKlogBmz5/Af5ROYpx1s3bi+12GZmfUVJ50OKUaF0TQ9USX93LnziV6GZGbWd5x0OqRAlREKAAQCYHh4Zy9DMjPrO046HVKICqNjSSf7OVIZ6WVIZmZ9x0mnQ4rUJR2ln046ZmZ7cNLpkGJUGE3JhpR8qlR6F5CZWR9y0umQ+paOBmrJx0nHzKyek06HFOuu6RQHSllh+FE4Zmb1nHQ6pEB1rHutUCwDILd0zMz24KTTIUUqVNLpHCzPBGDAScfMbA9OOh1SP5Bg5qz9AJC718zM9uCk0yFZSydLOvPmLgSyJ0+bmdluPU06kk6VdLektZIubLB+UNLVaf3Nkpam8hdJWiPp9vTzBXXbnJjK10r6qCTl8VmK7G7pHLDoEMBJx8xsvJ4lHUkF4DLgNOBo4CxJR4+rdg7wWEQcCXwEuCSVPwq8LCKeAZwNfL5um08AbwKWp9epXfsQdYqMjrV0Dn3KU4HsKQVmZrZbL1s6JwFrI+K+iBgGrgJOH1fndOCzafla4IWSFBH/GRG/TeV3AjNTq2gxMDciboqIAD4HvKLrnwQoRYVRZadz0UFLgKz1Y2Zmu/Uy6RwKPFj3fl0qa1gnIkaBrcD+4+q8CvhZRAyl+uta7BMASedKWi1p9caNGyf9IWrqr+kMFIsMRZGiWzpmZnuY1gMJJB1D1uX25r3dNiIuj4gVEbFi0aJFU46lSIXK2GNwYIgSxRid8n7NzH6X9DLprAcOq3u/JJU1rCOpCMwDNqX3S4CvAa+PiHvr6i9psc+uKFKhWjdmYZgSJScdM7M99DLp3AIsl7RMUhlYCawaV2cV2UABgDOBGyIiJM0HvgFcGBE/rlWOiA3ANkknp1Frrwe+3uXPAezZvQYwJCcdM7PxepZ00jWa84HrgV8A10TEnZIulvTyVO3TwP6S1gLvAGrDqs8HjgQuknRreh2Y1r0F+D/AWuBe4Jt5fJ76JxIAjFCkhJOOmVm9Yi8PHhHXAdeNK7uobnkX8OoG2/0N8DcT7HM1cGxnI21uxxM7GGSU6rhrOm7pmJntaVoPJOgXmx59hIKCSt01nSGVKDvpmJntwUmnAzZtzMYqVOuu6YxQohyeOdTMrJ6TTgds2bIJgGrd6RyWr+mYmY3npNMB27c+DkCVuiHTKjLolo6Z2R6cdDpg586UdFQ3ek0lym7pmJntwUmnA3YN7wQgxg2ZLjPcq5DMzPqSk04HDA8PAeOSjrvXzMyexEmnA0ZHshZNdY+kU3D3mpnZOE46HVCpZEkn6gYSjKrIIG7pmJnVc9LpgEq1Cjy5e63MCMPDvq5jZlbjpNMB1dqTB+qeSFChwKBGefi3uTzk2sxsWnDS6YAYm6xt9+kcTc9hW/fg/T2IyMysPznpdEI8uXut9vDPTY/8tuEmZmb7IiedDqgSAKju5tDa3DpbH9/Si5DMzPqSk05HZC2d+qRTe/jnzp3bexKRmVk/ctLpAKXutULdfDq1R+IMD+3sSUxmZv3ISacDVOteG6ibEy+yUzta9ZBpM7MaJ52OSC2dwu6kU1W2XK36BlEzsxonnQ5QSjrFQv3s39mprY4NpzYzMyedDhhI3WvF4uDuslpXW7reY2ZmTjodUWvpzJgxa6ysMHZ9xy0dM7MaJ50O2J10Zo6VFQvlPdaZmZmTTkcMpC60WbP3Gysrl2dk69zSMTMb46TTAbVrOnPnLRwrmzVrTrbO13TMzMY46XRArTUzf/4BY2XzUwIqeCI3M7MxTjodMBBZS+eAgw4ZK9v/gGy54Gs6ZmZjepp0JJ0q6W5JayVd2GD9oKSr0/qbJS1N5ftL+r6k7ZIuHbfND9I+b02vA7v9OWotnQULdx/qkKVPBaAQbumYmdUUW1fpDkkF4DLgRcA64BZJqyLirrpq5wCPRcSRklYClwCvBXYB7wWOTa/xXhcRq7v6AeoUqDISBUqF3c9eO+jgw9I6DyQwM6vpZUvnJGBtRNwXEcPAVcDp4+qcDnw2LV8LvFCSIuKJiPgRWfLpuUJUGaWwR5kKBXZFiZKfSGBmNqaXSedQ4MG69+tSWcM6ETEKbAX2b2Pfn0lda++V6uaQriPpXEmrJa3euHHj3kdfpxAVRsYlHYBhShTdvWZmNuZ3cSDB6yLiGcDz0uuPG1WKiMsjYkVErFi0aNGUDljgyS0dgCGKlJx0zMzG9DLprAcOq3u/JJU1rCOpCMwDNjXbaUSsTz8fB75I1o3XVQNRHZsptN6wyhTdvWZmNqaXSecWYLmkZZLKwEpg1bg6q4Cz0/KZwA0RaXxyA5KKkg5IyyXgpcAdHY98nAIVRtSoe61IyffpmJmN6dnotYgYlXQ+cD1QAK6IiDslXQysjohVwKeBz0taC2wmS0wASHoAmAuUJb0COAX4NXB9SjgF4LvAv3T7sxQbDCQAGKJEOTyfjplZTc+SDkBEXAdcN67sorrlXcCrJ9h26QS7PbFT8bWrQIXRBqdyWEXKvqZjZjbmd3EgQe6KUWnY0hl2S8fMbA9OOh2QtXSenHRG5Gs6Zmb1nHQ6oBjVhgMJhlRi0C0dM7MxTjodUGR0gpZOyS0dM7M6TjodMNE1nREKDOKWjplZTcukI+llkpycmihQpdrgFI2o6IEEZmZ12kkmrwXukfQhSU/vdkDTUbHJQAK3dMzMdmuZdCLij4DjgXuBKyX9JD0sc07Xo5smilFhtMFAglEnHTOzPbTVbRYR28imFrgKWAy8EviZpLd1MbZpo0iFSsOkU6DMCCMjTjxmZtDeNZ3TJX0N+AFQAk6KiNOAZwJ/0d3wpofSBN1roypQVoWHHxr/HFMzs31TO4/BOQP4SET8sL4wInZIOqc7YU0vRUYZbTCQoPbk6fW/uZclhy3NOSozs/7TTvfaQ+MTjqRLACLie12JapopRqXh1Aa1so0PPZR3SGZmfamdpPOiBmWndTqQ6axEhUqjlk66zvP444/lHZKZWV+asHtN0p8BbwGeKum2ulVzgB93O7DppEDjlk41le3Y+XjeIZmZ9aVm13S+CHwT+DvgwrryxyNic1ejmmaKVBreHFpNLZ2h4V15h2Rm1peaJZ2IiAckvXX8CkkLnXh2K1Gh0qCnMiIrq1SG8w7JzKwvtWrpvBRYAwSgunUBHNHFuKaViVo6pJZOper7dMzMoEnSiYiXpp/L8gtn+tmy+THmq/F9OrVxGhGVfIMyM+tT7dwcukrSWZJm5RHQdLNxY3bjZ6OWjgZSTnfSMTMD2hsy/WHgecAvJF0r6UxJM7oc17Tx2KaNAFQbnMrCQCktOemYmUEbTySIiBuBGyUVgBcAbwKuAOZ2ObZp4fEt2T04jZJOsVgGQFRzjcnMrF+18xgcJM0EXkY2zcEJwGe7GdR0sn37FqBx99pgaSbgpGNmVtMy6Ui6BjgJ+BZwKXBjRPhbNNm5cycAscfgvszMWfsB2c2jZmbWXkvn08BZ4SFYDQ0P7wAgGnSvzZu3EIABnzozM6D5Y3BeEBE3ALOB06U9/ycfEV/tcmzTwvDoEND4ms4BixYDbumYmdU0a+n8N+AGsms54wXgpANURkeBxi2dQw7P7p8tuKVjZgY0vzn0fWnx4oi4v36dJN8wmlRGs0fcNEo6By8+HMieWGBmZu3dp/OVBmXXduLgkk6VdLektZIubLB+UNLVaf3Nkpam8v0lfV/SdkmXjtvmREm3p20+qvH9gh1WSWMqGiUdFQrsihLFGO1mCGZm00azazpPB44B5kk6o27VXGDKN4em+34uI5uvZx1wi6RVEXFXXbVzgMci4khJK4FLyIZt7wLeCxybXvU+QXYv0c3AdcCpZE/L7opqtVL7QA3XD1Oi5O41MzOgeUvnKLIHfs4nu65Te51A9qU+VScBayPivogYBq4CTh9X53R23xN0LfBCSYqIJyLiR2TJZ4ykxcDciLgpIgL4HPCKDsQ6oUhdZ41aOgBDuKVjZlbT7JrO14GvS3p2RPykC8c+FHiw7v064FkT1YmIUUlbgf2BR5vsc924fR7aqKKkc4FzAQ4//PC9jX1MlttADe7TARhWiZKTjpkZ0N41nfMkza+9kbRA0hXdCykfEXF5RKyIiBWLFi2awp6yazpqNLUBWUunhJOOmRm0l3R+LyK21N5ExGPA8R049nrgsLr3S1JZwzqSisA8YFOLfS5psc8OqyWdRlMbpKTjlo6ZGdBe0hmQtKD2RtJC2nxmWwu3AMslLZNUBlYCq8bVWQWcnZbPBG6IWn9WAxGxAdgm6eQ0au31wNc7EOuEas9VGxhofEqGVaTspGNmBrSXPD4M/ETSl8lmDz0T+OBUD5yu0ZwPXA8UgCsi4k5JFwOrI2IV2SN4Pi9pLbCZLDEBIOkBspF0ZUmvAE5JI9/eAlwJzCQbtda1kWvpgwAwMEFLZ5gS5fDMoWZm0N7UBp+TtAb4g1R0xrhhzZMWEdeRDWuuL7uobnkX8OoJtl06QflqnjyMumtqLZ1SsfGpHFGRGTGcVzhmZn2trW6y1ALZSLo/R9LhEfGbrkY2TdSSTiHNnTPesIrMrT6RZ0hmZn2rnemqXy7pHuB+4EbgAbrdZTWNDKSkUy41TjojKlHG3WtmZtDeQIIPACcDv4qIZcALgZu6GtU0IrJrOuXyzIbrh1Vk0EnHzAxoL+mMRMQmslFsAxHxfWBFl+OaNgbSs9dmzpzTcP0oBY9eMzNL2rmms0XSfsAPgS9IegTwRYpEypLOfnPmNlw/oiKDeCCBmRm019I5HdgBXEA2ZfW9NJ5jZ59Ua+nMm3dAw/UjKlL2EwnMzID2hkzXWjVVdj9805Ja0llw4IEN11dUYJARRkZGKJVKeYZmZtZ32mnpWBOFNHrtwAMXN1w/QpGSKjyyYV3D9WZm+xInnSkaoEolxKz9Gl/TqaYnFTz463vzDMvMrC+1lXQkLZI0lUcx/84qRJWRJr2Uo2RJZ9PGh/MKycysb02YdJR5v6RHgbuBX0naKOmiibbZFxWojCWWRqppyoOt2zbnFZKZWd9q1tK5AHgO8F8iYmFELCCbZO05ki7IJbppoBDVpkmnklpBO3ZuzyskM7O+1Szp/DFwVkTcXyuIiPuAPyKbMsDIruk0Szq1aayHh3fkFZKZWd9qlnRKEfGkaaEjYiPgsb9JMZp3r0UaSDAy6htEzcyaJZ1m35L+Bk0KVBmdYC6dTNa9Fp7ewMys6c2hz5S0rUG5SFMcGBRatHSyWbYhqpW8QjIz61sTJp2IaPbfd0uKUWk6ZLpQyHoiJScdM7NJ3RwqyRO4Ja2GTJfLswBQOOmYmU32iQTqaBTTWDGaX9OZPWM/AAZw0jEzm2zSiY5GMY0VW7R05i/MHuRQCE/kZmY24cUISe+YaBWwX3fCmX5aJZ1DlxwBP8/qmZnt65qNXms8FWbmnzsdyHRViErT7rUjlh8DQMmzh5qZNR299td5BjJdFamwi/KE6/ebv4ChKFJ00jEza35NR9Jpkn4o6dH0ulHSS/IKbjootmjpAOxikLKv6ZiZNb2m8ybgzcC7gNWpeAXw95KWRMTlOcTX94pUqDS5pgOwixJlt3TMzJpe07kAeG5E1D+T/wZJpwE/Apx0SAMJWrR0dsotHTMzaN69pnEJB4CI2NSpg0s6VdLdktZKurDB+kFJV6f1N0taWrfu3an8bkkvrit/QNLtkm6VtHr8PjutnZbOECUGnXTMzJomnW2Snjm+MJU9PtUDSyoAlwGnAUcDZ0k6ely1c4DHIuJI4CPAJWnbo4GVwDHAqcDH0/5q/iAijouIFVONs5USFUbV/HanXRpk0A/8NDNr2r32F8AqSZ8B1qSyFcDZZHPqTNVJwNo0Rw+SrgJOB+6qq3M68P60fC1wqSSl8qsiYgi4X9LatL+fdCCuvVKM0bau6TjpmJk1aelExI/IZgodAN6QXgPAyWndVB0KPFj3fl0qa1gnIkaBrcD+LbYN4NuS1kg6d6KDSzpX0mpJqzdu3DjpD1GkQqVFS2dIZWbg7jUzs2YtHSLiIeCi2ntJBwAdu6bTJc+NiPWSDgS+I+mXEfHD8ZXS6LvLAVasWDHpx/oUqVBt8TShIZUYrLqlY2Y24belpJMl/UDSVyUdL+kO4A7gYUmnduDY64HD6t4vSWUN6yibmGYeWdKbcNuIqP18BPgaWbdb17TV0hkoMZOhboZhZjYtNPu2vBT4W+BLwA3An0bEwcDzgb/rwLFvAZZLWiapTDYwYNW4OqvIriEBnAncEBGRylem0W3LgOXATyXNljQHQNJs4BSyRNk17YxeG6HEDI9eMzNr2r1WjIhvA0i6OCJuAoiIX2bX8qcmIkYlnQ9cDxSAKyLiTkkXA6sjYhXwaeDzaaDAZrLERKp3Ddmgg1HgrRFRkXQQ8LUUXxH4YkR8a8rBTmDnzp0MttHSGVaJGW7pmJk1TTrVuuWd49Z1ZGqDiLgOuG5c2UV1y7uAV0+w7QeBD44ruw940jDvbtm6+VEOVlBtcXPoyECRGQwTEXQiYZuZTVfNks4zJW0jm8pgZlomvZ/R9cimgUc3buBgoBrNE8kwJcqqsO4397PkKUfkE5yZWR9q9pTp5v99N7ZufhSgZUunktbft/YuJx0z26dNduZQA554Imv8VVt0mY2m3L7xoXVdj8nMrJ856UzBE9uzpwFVo/lprLWEtm3f0u2QzMz6mpPOFOzctQOAaDF6raKspTO064mux2Rm1s+cdKZgeDgbBt3qiQTVdHlsZNRPJTCzfZuTzhSMVLKkEzS/pqPU0qn4oZ9mto9z0pmCykiWRFp1r2kgjW6LSrdDMjPra046U1CppCmoWwwkKAyUUz0nHTPbtznpTEE1JZFqi+61cim7l1Ya7XpMZmb9zElnCmKs5dL8NM6cNQeAgls6ZraPc9KZgkiPp1OLazoL5u0POOmYmTnpTEE2ywItH+K5eMlSAIq4e83M9m1OOlMRtQdxNz+NRx71ewAUw0nHzPZtTjpTsLt7rems38ydv5BdUaLspGNm+zgnnSlQSjoDA60fyD1EmaJnDzWzfZyTzhQMpKRTaCPp7KTMoJOOme3jnHSmQGkC1WKp1LLuLpUpuXvNzPZxTjpTUOteKxUGW9bdRZkZfvaame3jnHSmYCzplFvP3j2kMmV3r5nZPs5JZwoG0n06M2bMall3l9zSMTNz0pkCKWvpzN5vTsu6Q5QYxEnHzPZtTjpTMJBuDt1v9tyWdd3SMTNz0pmS2pDpOQv2b1l3eKDIjDZbOtetuoaf33rLlGIzM+tHzW+lt6ZqSWfRgYtb1h1WmZlttHS+/Y1reN6a/0l5zQh3/OsS1haXsH728bz1gr+ZcrxmZr3mls4UFFL32rx5bbR0KLZ1Teeuu37IHO3khvKJbC3M5oUjazhjyxe46nOXTTleM7Nec9KZggIVhqPArP32a1l3ZKDEzDaSzsLKI4zGAAf/4UU8573/jy8ufCPzeIKn3f8ZfnHbzzoRtplZz/Q06Ug6VdLdktZKurDB+kFJV6f1N0taWrfu3an8bkkvbnefnTQQVUbb7KEcUZGSKvx23f1N6x019Gt+MfAUjn/WcwF489v/lk/PPoMT4h5+teqv2LF9+5TjNjPrlZ5d05FUAC4DXgSsA26RtCoi7qqrdg7wWEQcKWklcAnwWklHAyuBY4BDgO9KelraptU+O6YQVUZp/dw1YCw5PXDPLzlkybKGdR584F5+r3ofqwafwzPqyt/2rk9y5d9v4Q27vsm1H30d5eWv4OWv/pOG+9ixfTs/uvEbbFh3P48/sYVqjCAKDKhAsVSmXJrBrJmzGZw5i8HBGRSKJQoDRSrVUSKCamWUaiWoVEbHpuNuZEDtfW4z62/HHvdsjlj+9NyO18uBBCcBayPiPgBJVwGnA/UJ4nTg/Wn5WuBSZTOmnQ5cFRFDwP2S1qb90cY+O6ZAhdE2G4uVdKoffnjdhHX+9auf5G0a5qHikwcmnHHex7nm42/kVUM/YOcdN3HlvdeztXAwil2UYxf7Vx5j2chveVr1QU7Rzsl9IDPb5/zT/a/nz9/1sdyO18ukcyjwYN37dcCzJqoTEaOStgL7p/Kbxm17aFputU8AJJ0LnAtw+OGHT+oDrC8/lS8VD+StbdStpJbBtm2bJqxTHn0EgMWLj3nSurnzF/Ka//WvfOrS97F0y495w65v7rF+e8zg7oHDuX7wWWwuzmeEMhXKhAaAKooqIhAVClRRVBAgZU9VIABBhMbektVoEGm08YnNbDpYuGBy33+Ttc8OmY6Iy4HLAVasWDGpb9G3v/MjbdetpG64XUM7JqyzbGgdD7KI1/zx+RPWefP5fw3AJ//5r9g5vI3B4jwWH7KUF774lZw4fyEnth2RmVn+epl01gOH1b1fksoa1VmnbHrOecCmFtu22mePZKd6eGSo4dod27dz3Ohabik9fY8PMJHz3u77dsxs+unl6LVbgOWSlkkqkw0MWDWuzirg7LR8JnBDREQqX5lGty0DlgM/bXOfPZF1c0G12njY9Je/eBkHagu/Lh/acL2Z2e+CnrV00jWa84HrgQJwRUTcKeliYHVErAI+DXw+DRTYTJZESPWuIRsgMAq8NSIbatVon3l/tkayhhpEtfGIsMe33JutLx2cW0xmZnnr6TWdiLgOuG5c2UV1y7uAV0+w7QeBD7azz34wMFDOftI46RwyuoFtMZNXvfYteYZlZpYrP5EgJ4Ol2uyijaesPmb4Pn5eOJKDD1mSX1BmZjlz0snJzBnZnDuFBi2dH3/v3ziKdfyqnO/QRTOzvDnp5GTuvIUAFBq0dG697T8AeHxgfp4hmZnlzkknJwcd8hQAig0eLVMZeRyAQrH1ZHBmZtOZk05Ojjwqe5paIZ7c0hmM7IbRw5cclWtMZmZ5c9LJyQGLDmRXlCg3SDrzqo+zIwZ50Wln9iAyM7P8OOnkaIgSpRh5UvnCylY2aCGz5szpQVRmZvlx0snRTgYpN0g6B1Ue46GB1rOPmplNd046OdqlcsPutcXVTTxSWNCDiMzM8uWkk6NdlBmMPZ+9tu4393KgtvBYwV1rZva7z0knR0MqP6l77cbrvwLAdu3Xi5DMzHLlpJOjXSozY1xLZ/OW3wIwqlm9CMnMLFdOOjkaosQgeyadgXSPzsxBX9Mxs999Tjo5atTSmRXbAXjm8c/tRUhmZrly0snRsErMGNfSmVd5nE0xh//6/D/sUVRmZvlx0snR8ECJmeNaOgdUtrBBvkfHzPYNTjo5Gqb4pGs6B1c281BhYY8iMjPLl5NOjkYGSsysSzoRwcGxiUcL83sXlJlZjpx0cjRCkZIqPLJhHQD/+ZMbmKudbPWNoWa2j3DSydGIigDcf+9dAPxszY0A7NDsnsVkZpYnJ50cVciSzkO//Q0AO4c2AxC+MdTM9hFOOjmqqADA9q2bAChVnwBg/tzFPYvJzCxPTjo5qpIlnSeGsmSzX3U7I1HgD174sl6GZWaWGyedHFVTSyeGtwKwoLKNh1nAsiM9TbWZ7RucdHI0e86R2dMHnvghwzue4MDKY2zw5G1mtg9x0snRueddwMdnvopnVO7jK5eew8HVzZ68zcz2KT1JOpIWSvqOpHvSz4bfvJLOTnXukXR2XfmJkm6XtFbSRyUplb9f0npJt6bXS/L6TO36H686n+9wIq/c8V0WxyY2Feb1OiQzs9z0qqVzIfC9iFgOfC+934OkhcD7gGcBJwHvq0tOnwDeBCxPr1PrNv1IRByXXtd18TNMylOXL+eWWb/PUJQoqsq2Ad8Yamb7jl4lndOBz6blzwKvaFDnxcB3ImJzRDwGfAc4VdJiYG5E3BQRAXxugu371v9610V8bOaZAGzX3B5HY2aWn14lnYMiYkNafgg4qEGdQ4EH696vS2WHpuXx5TXnS7pN0hUTddsBSDpX0mpJqzdu3DipDzEVr3zpm3nb4Ds55QWvy/3YZma9UuzWjiV9Fzi4war31L+JiJAUHTrsJ4APAJF+fhh4Y6OKEXE5cDnAihUrOnX8th1z7LF87Nhj8z6smVlPdS3pRMSEs5JJeljS4ojYkLrLHmlQbT3w+3XvlwA/SOVLxpWvT8d8uO4Y/wL8+2TjNzOzzutV99oqoDYa7Wzg6w3qXA+cImlB6iY7Bbg+dcttk3RyGrX2+tr2KYHVvBK4o1sfwMzM9l7XWjot/D1wjaRzgF8DrwGQtAI4LyL+NCI2S/oAcEva5uKI2JyW3wJcCcwEvpleAB+SdBxZ99oDwJu7/1HMzKxdygaA7dtWrFgRq1ev7nUYZmbTiqQ1EbFib7bxEwnMzCw3TjpmZpYbJx0zM8uNk46ZmeXGAwkASRvJRtFNxgHAox0Mp5P6NTbHtff6NbZ+jQv6N7Z+jQv2PranRMSivTmAk84USVq9t6M38tKvsTmuvdevsfVrXNC/sfVrXJBPbO5eMzOz3DjpmJlZbpx0pu7yXgfQRL/G5rj2Xr/G1q9xQf/G1q9xQQ6x+ZqOmZnlxi0dMzPLjZOOmZnlJyL8muQLOBW4G1gLXNilYxwGfB+4C7gTeHsqfz/ZPEK3ptdL6rZ5d4rpbuDFreIFlgE3p/KrgXKbsT0A3J6OvzqVLSSbWvye9HNBKhfw0XSM24AT6vZzdqp/D3B2XfmJaf9r07ZqM66j6s7LrcA24M97cc6AK8jmi7qjrqzr52iiY7QR2z8Av0zH/xowP5UvBXbWnbtPTjaGZp+zSVxd/90Bg+n92rR+aRtxXV0X0wPArXmfrxbfE33xt7ZHrN34otwXXkABuBc4AigDPweO7sJxFtf+IIA5wK+Ao9M/wnc2qH90imUw/eO6N8U6YbzANcDKtPxJ4M/ajO0B4IBxZR8i/QMHLgQuScsvIZuCQsDJwM11f7D3pZ8L0nLtH8ZPU12lbU+b5O/pIeApvThnwPOBE9jzi6rr52iiY7QR2ylAMS1fUhfb0vp64/azVzFM9DlbxNX13x3ZlCmfTMsrgatbxTVu/YeBi/I+Xy2+J/rib22PWPf2H7FfY7/kZ5NNKld7/27g3Tkc9+vAi5r8I9wjDrLJ8J49UbzpD+hRdn/R7FGvRSwP8OSkczewOC0vBu5Oy58CzhpfDzgL+FRd+adS2WLgl3Xle9Tbi/N1CvDjtNyTc8a4L6A8ztFEx2gV27h1rwS+0KzeZGKY6HO2OGdd/93Vtk3LxVRPzeKqKxfwILC8F+erwXFq3xN987dWe/mazuQdSvZHVrMulXWNpKXA8WRNf4DzJd0m6Yo0u2qzuCYq3x/YEhGj48rbEcC3Ja2RdG4qOyiy2V0ha2EcNMm4Dk3L48v31krgS3Xve33OIJ9zNNEx9sYb2T1BIsAySf8p6UZJz6uLeW9jmOy/nW7/7sa2Seu3pvrteB7wcETcU1fWk/M17nui7/7WnHSmCUn7AV8B/jwitgGfAJ4KHAdsIGva5+25EXECcBrwVknPr18Z2X99ogdxASCpDLwc+HIq6odztoc8ztFkjiHpPcAo8IVUtAE4PCKOB94BfFHS3G7GME7f/e7GOYs9/3PTk/PV4HtiyvvcG+0cw0ln8taTXbyrWZLKOk5SiewP6QsR8VWAiHg4IioRUQX+BTipRVwTlW8C5ksqjitvKSLWp5+PkF10Pgl4WNLiFPdisguvk4lrfVoeX743TgN+FhEPpzh7fs6SPM7RRMdoSdIbgJcCr0tfIkTEUERsSstryK6XPG2SMez1v52cfndj26T181L9plLdM8gGFdTizf18NfqemMQ+u/635qQzebcAyyUtS/+jXgms6vRBJAn4NPCLiPjHuvLFddVeCdyRllcBKyUNSloGLCe7ANgw3vSl8n3gzLT92WT9wa3imi1pTm2Z7NrJHen4ZzfY1yrg9cqcDGxNTfLrgVMkLUhdJqeQ9bFvALZJOjmdg9e3E9c4e/zvs9fnrE4e52iiYzQl6VTgXcDLI2JHXfkiSYW0fATZObpvkjFM9DmbxZXH764+3jOBG2pJt4U/JLveMdb9lPf5muh7YhL77P7fWrMLPn41f5GNAPkV2f9i3tOlYzyXrLl6G3XDRYHPkw1fvC390hfXbfOeFNPd1I34miheshE+PyUbCvllYLCNuI4gGxH0c7Ihmu9J5fsD3yMbPvldYGEqF3BZOvbtwIq6fb0xHXst8Cd15SvIvlzuBS6lzSHTadvZZP9LnVdXlvs5I0t6G4ARsn7wc/I4RxMdo43Y1pL16df+1mqjuV6Vfs+3Aj8DXjbZGJp9ziZxdf13B8xI79em9Ue0iiuVXwmcN65ubuerxfdEX/yt1b/8GBwzM8uNu9fMzCw3TjpmZpYbJx0zM8uNk46ZmeXGScfMzHLjpGNmZrlx0jGbJEn7S7o1vR6StD4tb5f08S4c70pJ90s6rwP7+ocU8zs7EZtZu4qtq5hZI5E95uQ4AEnvB7ZHxP/u8mH/MiKunepOIuIvJT3RiYDM9oZbOmYdJun3Jf17Wn6/pM9K+g9Jv5Z0hqQPSbpd0rfS87KQdKKypxGvkXT9uMe+THScKyV9QtJNku5Lx71C0i8kXZnqFFK9O9IxL+jqhzdrwUnHrPueCryA7InX/xf4fkQ8g2xmyf+eEs/HgDMj4kSyGSo/2Oa+F5DNCXMB2eNhPgIcAzxD0nFkLbFDI+LYdMzPdOpDmU2Gu9fMuu+bETEi6XayGS2/lcpvJ5vs6yjgWOA72bMUKZA946sd/xYRkfb9cETcDiDpzrTvG4EjJH0M+Abw7Y58IrNJctIx674hgIioShqJ3Q88rJL9GxRwZ0Q8e7L7Tvsaqiuvks2O+ZikZwIvBs4DXkP2QEeznnD3mlnv3Q0skvRsyOZFkXRMJ3Ys6QBgICK+AvwVcEIn9ms2WW7pmPVYRAxLOhP4qKR5ZP8u/4ns0fhTdSjwGUm1/2C+uwP7NJs0T21gNk2kEWn/3okh02l/7yefYd5mY9y9ZjZ9bAU+0KmbQ4E/AnyvjuXKLR0zM8uNWzpmZpYbJx0zM8uNk46ZmeXGScfMzHLz/wHYLkowgkLQCgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(model.BOLD.t_BOLD, model.BOLD.BOLD.T);\n",
    "plt.xlabel(\"Time [ms]\")\n",
    "plt.ylabel(\"BOLD activity\")"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}