{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# BE 240 Lecture : 2\n",
    "\n",
    "## Non-mass action propensities in bioscrape\n",
    "### _Ayush Pandey_\n",
    "\n",
    "### Biochemical Signalling Model\n",
    "\n",
    "A signal $S$ cooperatively binds to a transcription factor $F^0$ changing its conformation to $F^1$. $F^1$ binds to a gene $G^0$ activating it to $G^1$. A polymerase $P$ can bind to $G^1$ and then transcribes a transcript $T$. The transcript can bind to a ribosome $R$ and be translated to a protein $X$. The signalling molecule, transcripts, proteins are assumed to degrade via dilution. Cellularly machinery ($F$, $P$, $R$) and DNA ($G$) are assumed to be a constant concentration.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\\begin{align}\n",
    "2 S + F^0 &\\underset{k^u_1}{\\overset{k^b_1}\\rightleftharpoons} F^1\\\\\n",
    "F^1 + G^0 &\\underset{k^u_2}{\\overset{k^b_2}\\rightleftharpoons} G^1 \\\\\n",
    "G^1 + P &\\underset{k^u_3}{\\overset{k^b_3}\\rightleftharpoons} G^1 :P \\xrightarrow{k_{tx}} G^1 + P + T \\\\\n",
    "T + R &\\underset{k^u_4}{\\overset{k^b_4}\\rightleftharpoons} T :R \\xrightarrow{k_{tl}} T + R + X \\\\\n",
    "S &\\xrightarrow{\\delta} \\emptyset \\\\\n",
    "T &\\xrightarrow{\\delta} \\emptyset \\\\\n",
    "X &\\xrightarrow{\\delta} \\emptyset\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Using mass-action approximations, we can write a differential equation that describes the dynamics of each species. For example, if $C_1 = G^1:P$ and $C_2 = T:R$, we have the following mRNA dynamics using mass-action propensities:\n",
    "\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{dT}{dt} = k_{tx}C_1 + (k_{tl} + k_4^u) C_2 - k_4^b T R - \\delta T\\\\ \n",
    "\\\\\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "* **How do we model the fact that some binding/unbinding reactions are fast reactions compared to the catalysis reactions?**\n",
    "\n",
    "\n",
    "* **Moreover, if we are only interested in the dynamics of the protein production, can we write simpler models that express the effective dynamics appropriately?**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Reduced model representations\n",
    "\n",
    "## Time-Scale Separation\n",
    "\n",
    "* **QSSA : Quasi steady-state approximation**. In a dynamic model, if a species exhibits a different intrinsic time scale and tends to reach an equilibrium state quicker than others, then we approximate such a species by its quasi steady state. For example, in the above example, if the binding and unbinding reaction of the transcript with the RNAP are fast, we can assume that the complex $C_1$ is at QSS. We can then write,\n",
    "\\begin{align}\n",
    "\\frac{dC_1}{dt} = 0\n",
    "\\end{align}\n",
    "On solving the above, we get an algebraic relationship for $C_1$ as function of other species that we can substitute into all other equations to eliminate $C_1$ dynamics from the model.\n",
    "* Formally, singular perturbation theory is a mathematical tool to derive QSSA conditions and steady state expressions. For a rigorous mathematical treatment of time-scale separation as applied to biological systems, refer to Chapter 3 of [[1]](http://www.cds.caltech.edu/~murray/books/AM08/pdf/bfs-dynamics_14Sep14.pdf) and for more information on singular perturbation theory refer to [[2]](https://books.google.com/books/about/Nonlinear_Systems.html?id=t_d1QgAACAAJ)\n",
    "\n",
    "[1]: Del Vecchio, Domitilla, and Richard M. Murray. Biomolecular Feedback Systems. Princeton, NJ: Princeton University Press, 2015.\n",
    "\n",
    "[2]: Khalil, Hassan K., and Jessy W. Grizzle. Nonlinear Systems. Vol. 3. Upper Saddle River, NJ: Prentice hall, 2002."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Conservation Laws\n",
    "\n",
    "If a species in a model is conserved in amount, but undergoes chemical reactions, we can write a conservation law that relates the total amount of that species with the various amounts in which that species is present in the system. We can use these conservation laws to eliminate redundant dynamics from our system models. We will demonstrate these techniques using a simple example."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## A Simple Example: An Enzymatic Reaction System\n",
    "\\begin{align}\n",
    "    E + S &\\underset{d}{\\overset{a}\\rightleftharpoons} C \\xrightarrow{k} P + E \n",
    "\\end{align}\n",
    "Here, $C$ is the complex formed when the enzyme (E) binds to the substrate (S). P is the product resulting from the modification of the substrate S due to the binding with the enzyme E. Refer to [1] for more information.\n",
    "The full mass-action kinetics based model with four species is given by:\n",
    "\\begin{align*}\n",
    "\\frac{dS}{dt} &= -aES + dC \\\\\n",
    "\\frac{dC}{dt} &= aES - (d+k) C\\\\\n",
    "\\frac{dE}{dt} &= -a ES + dC + kC \\\\\n",
    "\\frac{dP}{dt} &= kC\n",
    "\\end{align*}\n",
    "\n",
    "The total enzyme concentration is usually assumed to be constant. Hence, we get the following **conservation law**:\n",
    "\\begin{align}\n",
    "E + C = E_{tot}\n",
    "\\end{align}\n",
    "We can eliminate E from the system dynamics by substituting $E = E_{tot} - C$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**Time-Scale Separation Assumption**:\n",
    "\n",
    "\n",
    "Next, we assume that the binding and unbinding reactions are at a faster time-scale than the catalysis reaction. Hence, we assume QSSA for C:\n",
    "\\begin{align}\n",
    "\\frac{dC}{dt} &= 0 \\Rightarrow a(E_{tot} − C)S − (d+k)C = 0\\\\\n",
    "\\Rightarrow C &= \\frac{E_{tot} S }{S+ K_m}, \\quad K_m = \\frac{d + k}{a}\\\\\n",
    "\\\\\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "The reduced model then becomes:\n",
    "\\begin{align}\n",
    "\\frac{dP}{dt} &= rate . \\frac{S}{S+ K_m}\\\\\n",
    "\\\\\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "This is a positive hill function with $n = 1$, $K = K_m$, and rate$ = k E_{tot}$.\n",
    "The mathematical conditions on the parameters under which this assumption is justified are discussed in [1]. \n",
    "\n",
    "\n",
    "[1]: Del Vecchio, Domitilla, and Richard M. Murray. Biomolecular Feedback Systems. Princeton, NJ: Princeton University Press, 2015"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "**All other kinds of Hill functions can be derived in a similar fashion. Near the end of BE240, we will discuss tools to do these model reduction techniques using an automated package.**\n",
    "\n",
    "\n",
    "\n",
    "**Here is a list of the Hill function propensities available in bioscrape:**\n",
    "\n",
    "1. Positive Hill \n",
    "\\begin{align}\n",
    "rate \\times \\frac{X^n}{K^n + X^n}\\\\\n",
    "\\\\\n",
    "\\end{align}\n",
    "\n",
    "1. Negative Hill \n",
    "\\begin{align}\n",
    "rate \\times \\frac{K^n}{K^n + X^n}\\\\\n",
    "\\\\\n",
    "\\end{align}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "\n",
    "1. Positive Proportional Hill \n",
    "\\begin{align}\n",
    "rate \\times d \\times \\frac{X^n}{K^n + X^n}\\\\\n",
    "\\\\\n",
    "\\end{align}\n",
    "\n",
    "1. Negative Proportional Hill \n",
    "\\begin{align}\n",
    "rate \\times d \\times \\frac{K^n}{K^n + X^n}\\\\\n",
    "\\\\\n",
    "\\end{align}\n",
    "\n",
    "Here $X$ and $d$ are species identifiers and $K$ and $rate$ are parameters. Detailed description of non-massaction propensities is available on the [bioscrape wiki](https://github.com/ananswam/bioscrape/wiki/Propensities).\n",
    "\n",
    "An [interactive tool](https://www.physiologyweb.com/calculators/hill_equation_interactive_graph.html) online to practice your Hill functions!\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Assignment: Try to write bioscrape models for the enzymatic reaction system - \n",
    "* full dynamics\n",
    "* reduced dynamics with a Hill function. \n",
    "\n",
    "### Compare the results and find the parameter regime in which the reduced model works. Solution available in Chapter 3 of [1].\n",
    "\n",
    "\n",
    "\n",
    "[1]: Del Vecchio, Domitilla, and Richard M. Murray. Biomolecular Feedback Systems. Princeton, NJ: Princeton University Press, 2015"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "## Example 2 : Using bioscrape non-mass action propensities \n",
    "Consider the transcription-translation model where $G$ is a gene, $T$ is a transcript, $X$ is a protein.\n",
    "We can write the following reduced order dynamics:\n",
    "1. $G \\xrightarrow[]{\\rho_{tx}(G, S)} G + T$; \n",
    "\\begin{align} \n",
    "\\rho_{tx}(G, S) = G k_{tx}\\frac{S^{2}}{K_{S}^{2}+S^{2}}\n",
    "\\\\\n",
    "\\end{align}\n",
    "Here, $S$ is the inducer signal that cooperatively activates the transcription of the gene $G$. Since, this is a positive activation of the gene by the inducer, we have a positive proportional Hill function.\n",
    "\n",
    "2. $T \\xrightarrow[]{\\rho_{tl}(T)} T+X$; \n",
    "\\begin{align} \n",
    "\\rho_{tl}(T) = k_{tl} \\frac{T}{K_{R} + T}\n",
    "\\\\\n",
    "\\end{align}\n",
    "The parameters $k_{tl}$ and $K_R$ model effects due to ribosome saturation. This is the similar Hill function as derived in the enzymatic reaction example. \n",
    "\n",
    "3. $T \\xrightarrow[]{\\delta} \\varnothing$; massaction kinetics at rate $\\delta$.\n",
    "4. $X \\xrightarrow[]{\\delta} \\varnothing$; massaction kinetics at rate $\\delta$.\n",
    "\n",
    "The first reaction uses a proportional positive hill function as its rate function to represent induction. The second reaction uses a positive hill function function to represent ribosome saturation. The third and fourth reactions reaction represent degredation via dilution. This model is constructed below and simulated both stochastically and deterministically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "from bioscrape.simulator import py_simulate_model\n",
    "from bioscrape.types import Model\n",
    "\n",
    "#Create a list of species names (strings)\n",
    "species = [\"G\", \"T\", \"X\", \"S\"]\n",
    "\n",
    "#create a list of parameters in the form (param_name[string], param_val[number])\n",
    "params = [(\"ktx\", 2), (\"ktl\", 5), (\"KI\", 10), (\"n\", 2.0), (\"KR\", 20), (\"delta\", .1)]\n",
    "\n",
    "#create reaction tuples in the form:\n",
    "#(Inputs[string list], Outputs[string list], propensity_type[string], propensity_dict {propensity_param:model_param})\n",
    "rxn1 = ([\"G\"], [\"G\", \"T\"], \"proportionalhillpositive\", {\"d\":\"G\", \"s1\":\"S\", \"k\":\"ktx\", \"K\":\"KI\", \"n\":\"n\"})\n",
    "rxn2 = ([\"T\"], [\"T\", \"X\"], \"hillpositive\", {\"s1\":\"T\", \"k\":\"ktl\", \"K\":\"KR\", \"n\":1}) \n",
    "\n",
    "# Or alternatively, use general propensity (simulation will be slower):\n",
    "# rxn2 = ([\"T\"], [\"T\", \"X\"], \"general\", {\"rate\":\"ktl * T/(KR + T)\"})\n",
    "\n",
    "#Notice that parameters can also take numerical values instead of being named directly\n",
    "rxn3 = ([\"T\"], [], \"massaction\", {\"k\":\"delta\"})\n",
    "rxn4 = ([\"X\"], [], \"massaction\", {\"k\":\"delta\"})\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "C:\\Users\\apand\\Anaconda3\\lib\\site-packages\\ipykernel_launcher.py:8: UserWarning: The following species are uninitialized and their value has been defaulted to 0: T, X, \n",
      "  \n"
     ]
    }
   ],
   "source": [
    "#Create a list of all reactions\n",
    "rxns = [rxn1, rxn2, rxn3, rxn4]\n",
    "\n",
    "#create an initial condition dictionary species not included in the dictionary will default to 0\n",
    "x0 = {\"G\":1, \"S\":200}\n",
    "\n",
    "#Instaniate the Model object\n",
    "M = Model(species = species, parameters = params, reactions = rxns, initial_condition_dict = x0)\n",
    "\n",
    "#Simulate the Model deterministically\n",
    "timepoints = np.arange(0, 150, .1)\n",
    "results_det = py_simulate_model(timepoints, Model = M) #Returns a Pandas DataFrame\n",
    "\n",
    "#Simulate the Model Stochastically\n",
    "results_stoch = py_simulate_model(timepoints, Model = M, stochastic = True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAr8AAAEWCAYAAABv4v9VAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOydd3gbRdrAfyNZbrETp/dKem9ASCgJ5LjQA6HDUY56HHBw9DvuLhxwcJSDL8ARQi85egudEFJJSEjvnfRe7NiOm6T5/phda7VeWZIt2bI1v+fRoy2zu6Pd0cy777xFSCnRaDQajUaj0WiSAVdtV0Cj0Wg0Go1Go6kptPCr0Wg0Go1Go0katPCr0Wg0Go1Go0katPCr0Wg0Go1Go0katPCr0Wg0Go1Go0katPCr0Wg0Go1Go0katPCrSTqEEH8TQkys7XpoNBqNJhghxBVCiO9rux6a+o3QcX41VUUIUWBZzQRKAJ+xfpOUcnLN1yq2CCEeAdpJKa9x2OcG8iyb7PfgOinl+3GvpEaj0cQQIcQWoCWqLysEvgZuk1IWVHZciHONB7pKKa+MZR1DXOtpoK+U8reWbc8a1z873tfX1B205ldTZaSUWeYH2AacY9lWQfAVQqTUfC3jVwcppc92D3YBZ1i2acFXo9HUVc4x+rXBwLHAg/YCQpFIcsTfgGOEENcCCCFOAK4Gbq7VWmkSjkRqtJp6hhDiESHE+0KId4UQ+cCVQogThBA/CyFyhRC7hRAThBAeo3yKEEIKIW4SQmwUQhwWQkywnK+7EGKWECJPCHFACPE/y75+QogfhBCHhBB7hBD3VlKHR4QQbxj7uxrXvEEIscv43GnsOxu4F7hCCFEghFhUYzdPo9FoEgAp5U7gG6AvgBBihhDiUSHET8BRoIsQoo0QYorR/24UQtxglB0D/AW4xOhDlxnbGwkhXjXGgJ1Gn+w29l0jhJhjXt/on28WQmwwxoQXhBAiRF2PAtcDTwkhOgGvAfdLKXfE5+5o6iq1ronT1HvOB8YBVwBpqA70T8AioAPwLbAeeN5yzJnAEKAxsFgIMUVK+QPwKPAVMBJINcoghGgE/AD8CzjLuE7PSurwV4d6ngx0BboBPwohlkgpvxRCPEEIsweNRqOp7wgh2qP65E8sm38HnAGsAwSq/10FtEH1vVOFEJullN8KIf5FRbOHN4G9qD63AfAlsB14KUQ1zkZpnxuixo4vUGNHBaSUM4QQHwELgZXApGh/s6b+ozW/mngzR0r5hZTSL6UsklL+IqWcL6X0Sik3ozqmU2zHPCalzJNSbgFmAAON7WVAJ6C1lLJYSvmTsf1cYLuU8v+klCVSyiNSygWh6hCing9JKY9KKZehOubLqv3LNRqNpu7ymRAiF5gDzEQpF0zekFKuklJ6gVbAicB9Rr+8FHgFJSBXQAjREiU43yGlLJRS7gOeAS6tpC6PSylzpZTbgOkExoRQzAaaApOldmzSOKCFX0282W5dEUL0FEJ8ZZgmHAH+CTSzHbPHsnwUyDKW7wI8wEIhxAohxNXG9vbAxkjrEEGZrSgNhkaj0SQrY6WUOVLKjlLKW2yKA2t/2QY4JKXMt2zbCrQNcd6OqH58t2H+lovS+LaopC6hxoQKCCGaAk8BzwL/FELkVHJeTZKihV9NvLG/db+EmorqKqVsCPwdNW0W/kRS7pZSXi+lbA38EZgkhOiM6oiPiaIOTrS3LHdAOa9FeqxGo9EkE9Z+cRfQRAiRbdnWAdjpUBZUf10CNDOE6xwpZUMpZZ8Y1e1Z4Fsp5Z3ALJQgrNEEoYVfTU2TjQoPViiE6AXcFOmBQoiLhRCmNiEX1an6gClAByHErUKIVCFEQyHEcVHW629CiAwhRD+Ud7AZqWEv0CmUg4VGo9EkM1LK7cBc4DEhRLoQoj9wHWBG/DH7UJdRfjfwPfC00Ve7hBDHCCHs5m9RI4Q4E/gN8Gdj023AWCHEqOqeW1O/0MKvpqa5CyVc5qO0wNGEAzse+EUIUYhyvvijlHKblDIP1eGNA/ahHOii7UjnAJtRnfJjUsofje3vo5zrDgkhFoQ6WKPRaJKYy1D+GLuAT4F/SCmnGvs+NL4PCiEWG8tXofrV1cBh4COgdXUqYGieJwK3SykPARj2xHcBLwshMqpzfk39Qie50CQ1QoiuwAYppdbsajQajUaTBGjNr0aj0Wg0Go0madDCr0aj0Wg0Go0madBmDxqNRqPRaDSapEFrfjUajUaj0Wg0SUONpjdu1qyZ7NSpU01eUqPRaGLCokWLDkgpm9d2PWoS3WdrNJq6TKh+u0aF306dOrFw4cKavKRGo9HEBCHE1tquQ02j+2yNRlOXCdVva7MHjUaj0Wg0Gk3SoIVfjUaj0Wg0Gk3SoIVfjUaj0Wg0Gk3SUKM2vxqNJpiysjJ27NhBcXFxbVdFY5Cenk67du3weDy1XZWERLfZxEO3WY0mOrTwq9HUIjt27CA7O5tOnTohhM6wXNtIKTl48CA7duygc+fOtV2dhES32cRCt1mNJnrCmj0IIdKFEAuEEMuEEKuEEA8Z298QQvwqhFhqfAbGv7oaTf2iuLiYpk2baiEiQRBC0LRp0zqt1Yx3n63bbGJRH9qsRlPTRKL5LQFOlVIWCCE8wBwhxDfGvnuklB/Fr3oaTf1HCxGJRT14HnHvs+vBPapX6Oeh0URHWM2vVBQYqx7jo3Mia6qGtwSWvAM6rbZGExd0n63R1GNWfQqz/wOlR2u7JnWaiKI9CCHcQoilwD5gqpRyvrHrUSHEciHEM0KItBDH3iiEWCiEWLh///4YVVtTZ5n5BHz+R1j7ZW3XRBOC8ePH89RTT4Xc/9lnn7F69eq4Xf/6668Pe/6JEyfy1ltvhdw/Y8YM5s6dG3H5+kay9dm6zWqShg+vgWkPwZY5tV2TOk1EDm9SSh8wUAiRA3wqhOgLPADsAVKBScB9wD8djp1k7Gfo0KFa+5DsFOxR38V5tVsPTZX57LPPOPvss+ndu3fEx3i9XlJSIvOvfeWVV8KWufnmmyvdP2PGDLKyshg+fHhE5esbus8ORrdZTb3AOmNapjW/1SGqOL9SylxgBjBGSrnbmF4rAV4HjotD/TT1DfO/q80eEopHH32UHj16MHr0aNatWwfApk2bGDNmDEOGDOGkk05i7dq1zJ07lylTpnDPPfcwcOBANm3a5FgO4JprruHPf/4zo0aN4r777mP8+PFcffXVnH766XTq1IlPPvmEe++9l379+jFmzBjKysoAGDlyZHlK3aysLP76178yYMAAhg0bxt69e4FgTd+ECRPo3bs3/fv359JLL2XLli1MnDiRZ555hoEDBzJ79uyg8hs3bmT06NEMGDCAwYMHs2nTphq91zVJfe6zdZvVJB2+0sCyVzs4Voewr7VCiOZAmZQyVwiRAYwG/i2EaC2l3C2Upf1YYGWc66rR1Gse+mIVq3cdiek5e7dpyD/O6VNpmUWLFvHee++xZMkSvF4vgwcPZsiQIdx4441MnDiRbt26MX/+fG655RZ+/PFHzj33XM4++2wuvPBCAE477TTHcgDr16/nhx9+wO12M378eDZt2sT06dNZvXo1J5xwAh9//DFPPPEE559/Pl999RVjx44NqlthYSHDhg3j0Ucf5d577+Xll1/mwQcfDCrz+OOP8+uvv5KWlkZubi45OTncfPPNZGVlcffddwMwbdq08vJXXHEF999/P+effz7FxcX4/f5q3+dEoib7bN1mdZvV1CBWgVcLv9Uikjmd1sCbQgg3SlP8gZTySyHEj0YnK4ClgJ6j0WjqILNnz+b8888nMzMTgHPPPZfi4mLmzp3LRRddVF6upKSkwrEFBQWVlrvoootwu93l62eccQYej4d+/frh8/kYM2YMAP369WPLli0Vzp+amsrZZ58NwJAhQ5g6dWqFMv379+eKK65g7NixFQQRO/n5+ezcuZPzzz8fUMkB6iH1vs/WbVaTlJQVOy9roias8CulXA4Mcth+alxqpKn7rP8efvgHXDAJWvULbH//SljzhVqecis0bANdT6udOiYg4bRd8cQeKsnv95OTk8PSpUsrPS5cuQYNGgStp6UpHyuXy4XH4ym/rsvlwuv1VjjeWsbtdjuW+eqrr5g1axZTpkzh4YcfZtWqVSHrK5PA3KYm+2zdZnWbjSsLX4Nt8+GCl2q7JomB1vzGjKhsfjWaiFj+PuxbDbssg4uUAcHXZPXnNVsvjSMnn3wyn376KUVFReTn5/PFF1+QmZlJ586d+fDDDwE1AC9btgyA7Oxs8vPzAWjYsGHIcjWB3+9n+/btjBo1iieeeILc3FwKCgqC6milYcOGtGvXjs8++wxQGr+jR7XjSF1Dt9kkabNf3gnL36vtWiQO3hLnZU3UaOFXE3tMjYz1zdRXFrqcplYZPHgwl1xyCQMHDmTcuHGcdNJJAEyePJlXX32VAQMG0KdPHz7/XL2sXHrppTz55JMMGjSITZs2hSxXE/h8Pq688kr69evHoEGDuPPOO8nJyeGcc87h008/LXcesvL2228zYcIE+vfvz/Dhw9mzZ0+N1VcTG3Sb1W02KfEWOS9rokbU5JTK0KFDpekRq6nHfHQdrPwITn8Eht+mthXnweMdgssNuQbO+b8ar14isWbNGnr16lXb1dDYcHouQohFUsqhtVSlWsGpz9ZtNjGpl89lfCPjW4fGBGD7Anj1N2p52C0w5rHarU8dIFS/rTW/mvgRZJ/kMEUjtceyRqPRaMKgo1sotM1vzNDCryb2SJ/6zt0GJQXqc/RgxXK7l8HORWo/gM8LhzbXXD01Go1Gk/gc2VHbNUgMrBEedi+HsiLlW1N0uPbqVEfRwq8m9pha3sVvKVOHx9rC9w9WLLd7Gbx8Knxyg1qf+W+YMAj2r6u5umo0Go0msXm2X/gyyYBV27tzITzaCiadAv+7tPbqVEfRwq8m9pQZhvgp6QEtsN8HrhS4+ku48pPg8od+Vd+bp6vvwgM1U0+NRqPRJCZ+n21dmz6EjPBwZFfN1qMeEFnico0mGsw/qN8S39LvhXbHQeeTKpZ32ZqhT4dw0Wg0mqTGHiHI7wVXau3UJVEIFeHB7xBNSVMpWvOriT3m1EyQ8OsDl9u5vNsm/Or4hRqNRpPc2AU6f8VkIUlHqLHRV1qz9agHaOFXE3ucvFD93ooaXhP79jIdv7C2efbZZ6scSH/8+PE89dRT1a7DG2+8wa5dgem866+/ntWrV1f7vJr6iW6z9YwKml+t3Qw5Nvr0i0G0aOFXEzv8Ppjxbziys+K+A+tDC787foEPrlbfAMveg/Xfxa+emrBUR5CIFXZB4pVXXqF37961WCNNIqPbbAJzZDcsmRzdMRWEX59zuWQilOa3JA9+ne28T+OIFn41sWP/OpjxL9VpteofvK84N1j4Hf1Q8P7VnwWWN3wH/7s4fvXUBFFYWMhZZ53FgAED6Nu3Lw899BC7du1i1KhRjBo1CoB3332Xfv360bdvX+67777yY7/99lsGDx7MgAEDOO2008q3r169mpEjR9KlSxcmTJhQvn3s2LEMGTKEPn36MGnSJEBlvLrmmmvo27cv/fr145lnnuGjjz5i4cKFXHHFFQwcOJCioiJGjhyJmXAh1HU1yYFus3WMyRfB57fA0UORH2OfytdmD2pWVbjgpLsr7pvzn5qvTx1GO7xpYofpqDbuVWjcEV4cHrzfavN74h3qY2bw0cA398OeFbE9Z6t+cMbjlRb59ttvadOmDV999RUAeXl5vP7660yfPp1mzZqxa9cu7rvvPhYtWkTjxo05/fTT+eyzzxgxYgQ33HADs2bNonPnzhw6FBjY1q5dy/Tp08nPz6dHjx784Q9/wOPx8Nprr9GkSROKioo49thjGTduHFu2bGHnzp2sXLkSgNzcXHJycnj++ed56qmnGDo0ODnP/v37Q15XU8PoNgvoNhsWczYwGgHWNHPoMBy2zdXCLyjhNyUDTvub+pjjZ8cRFTXlmkrRml9N7DD/fO5UFebMTiizB02t0q9fP3744Qfuu+8+Zs+eTaNGwS8kv/zyCyNHjqR58+akpKRwxRVXMGvWLH7++WdOPvlkOnfuDECTJk3KjznrrLNIS0ujWbNmtGjRgr179wIwYcIEBgwYwLBhw9i+fTsbNmygS5cubN68mdtuu41vv/2Whg0bVlrfyq6rSQ50m62jROPMbI4nnvTg9WTGWwwpaRW3u1K001uUaGlEEzvKhV+Ps/Dr9tRsfeoaYbRd8aJ79+4sWrSIr7/+mgceeIDTTz89aL+U0vE4KSVCCMd9aWmBDtrtduP1epkxYwY//PAD8+bNIzMzk5EjR1JcXEzjxo1ZtmwZ3333HS+88AIffPABr732Wsj6VnZdTQ2j26xus5Fg/vaohF9DmPNkqm+t+VXCryej4nZ3KpTk13x96jBa86uJHWZn5fY4/0Gj1fxqB4caYdeuXWRmZnLllVdy9913s3jxYrKzs8nPV53p8ccfz8yZMzlw4AA+n493332XU045hRNOOIGZM2fy668qSUm4qdy8vDwaN25MZmYma9eu5eeffwbgwIED+P1+xo0bx8MPP8zixYsBgupgJdrrauofus3WMcyXkVBxap0o1/waY4keD1R6YyfNr9ujNeNRojW/mtjw+a2w7mu1HCuzh9yt0KRL9eumqZQVK1Zwzz334HK58Hg8vPjii8ybN48zzjiD1q1bM336dB577DFGjRqFlJIzzzyT8847D4BJkyZxwQUX4Pf7adGiBVOnTg15nTFjxjBx4kT69+9Pjx49GDZsGAA7d+7k2muvxW9kcHrssccAuOaaa7j55pvJyMhg3rx55edp3rx5VNfV1D90m60jrPtGjQ1FhrD/4yNw+fuRHVsu/Bqa3xeOhbZD4bqp4HLQ2017GFZ/DjdOh7Ts6tc90cjbASs/guY9K+5ze7TZQ5SIUNND8WDo0KHS9HzV1DOsjms3zoQ2A2Hu89CsG8x5Flr0gkFXQtvBwcftXg6znoDNs1S4lqxWUHIEyo7CNV9BpxNr9nfUMGvWrKFXr161XQ2NDafnIoRYJKUcGuKQeolTn63bbGKSkM/lx0dg1pOB9dYD4aaZkR27ZQ68cRZc/BZ8cFVg+192QWqDiuXNMejWhWrcqW9smg5vj4UTboXfPqq2mb+530WwcxHcvqT26peghOq3teZXE3vcRgrK4beq7+6/DV22dX+45B347BZYOhlOvhta9oHXz9DTOBqNRlOXsSdliMZu19RkNmhhO2exs/Bbflw9HTdMe+m+4yruc2mzh2gJa/MrhEgXQiwQQiwTQqwSQjxkbO8shJgvhNgghHhfCJHkSbc15VTHsS0lLWAeoW28NJqo0X22JmGwO7hFMzVvZi1z25qpUwZRK/XVMc783aGcybXwGxWROLyVAKdKKQcAA4ExQohhwL+BZ6SU3YDDwHXxq6amTlEd4Ve4LMJvPe3EbNSk6ZEmPPXgecS9z64H96hekbDPwy6oRiOglTtQ2yaok1349YQSfrXNbzSEFX6losBY9RgfCZwKfGRsfxMYG5caahKbosOwf33wNlc1hF8pk0r4TU9P5+DBg4k7eCUZUkoOHjxIerrDAFNHiHefrdtsYpGwbTZ/L+RuC97mK4NdS8AbQlArOgwHNqhlM3mKXfO7/lvI3R68zdoWQ80Y+spg2/zQ1050KtX8pqrfd2gzbPoRinLV9vy9UHyk5upYh4jI5lcI4QYWAV2BF4BNQK6U0pROdgBtQxx7I3AjQIcOHapbX02i8e7lKvuOFXtnFQmtByib34Ztkkr4bdeuHTt27GD//v21XRWNQXp6Ou3atavtalSLePbZus0mHgnZZj+5Hn6dFbztyA6YNBJG/gVG3lfxGHM8uXMVzDRiSKfnQPNesH+NWv/uL7D4bfjjz4HjtgUia5RnhrOz5B348g747b/ghD9W+WfVGmUOwm/TrnBwo9L8+stgwqDAvvF58HR3yOkAd8Q4C2M9ICLhV0rpAwYKIXKATwEnl1JHNYCUchIwCZTncBXrqUlU7IIvVJymioRjr4fOJ6swLuabfxIIvx6Ppzzrk0YTK+LZZ+s2q4mI4jxodxzsWKDWe54Na79Uy7tCRCUwx5O8nYFtDVvDtV+rY6fcprYV2eIkH7Wshxo3Cg9ULFuXcNL83jRbbZ/3QmizB7v2XQNEmeRCSpkLzACGATlCCFPKaQfsim3VNHWWqmh+XW4VDk0ItQxJIfxqNPFE99maWsPvgyxLpIbMKFI625NhZDZRGkyTMpvdr9UOONS4YWqEo403nyiYzoPWJBepmereuFNB+munXnWUSKI9NDe0BwghMoDRwBpgOnChUexq4PN4VVJTx6iK8GslicweNJpYo/tsTULg9wYUGQBpDSM/1h4iDYKFVrvTm3XdF2LcqOvRELxFyp/Gek9NqjLbmuREcsdaA28aNmQu4AMp5ZdCiNXAe0KIR4AlwKtxrKemLlHdN2szWoQWfjWaqqD7bE3t4/cGjwVOKe9DUZxXcZv1XL4S5eQmhFqPRPNrmgXU1XHFWxL6HjopnHSo0EoJK6VIKZcDgxy2bwaOi0elNHUcs0OqKmYn9+0D0H0M/PwitD8eep5Z/bppNPUc3WdrEgJfWbDA6hQFaN03MPtpaDNYZQU1+eyWimWFTeO57D0YeJlatppBbJ5ecawo2AfznlfLOxbA1/fAyAeiM8WobbbMCa1YchJ+dy8LLK/7BnqcEZ96VYeyIuWkd9wNcNJdNXrpqGx+NZogrCFUmnVXHVj3GPzBMhpDo/YqxfHit+GnZ+G9y6p/Xo1Go9HUDH6fEtbOfQ66jg6O/25qMN+9FHb8Agtegs/+ENgvDa3lwCtCn3/+i4Flq+b34KaKZTdbUir/OgsWTAqOEFEXkBJKC533tRkETbqoyBgms58OLM97Ib51qypzn4f83TDtnzV+aW0ooqk6Zodz1tMqWkOscHvg99/CM32gND9259VoNBpNzWCaPQy+Sn3mPhfYJ0NMyXc8EbbOUcvjHUwfANoOgaxWwVEMvCWAgA7DnKMe2B3oIHyyjETDVxJ69rP9cXC7EUGjtBD+1SZgOpLeKHF/a2lBYNlXVr0EWVGiNb+aqmM6JTgF3a4u5jlLCiovp9FoNJrEw27za52at0drMElvFNm5Xe5g211vkdImp6Q7C3pO1wtVh0SlrDiysdZtRIMwE12kN6obv7WGBXQt/GqqTnnolTgKv6Va+NVoNJo6RwXh16LVCyXoZOQ4b4dgXxJXSnAyC2+JCgEWSvgtTw3cIHwdEhVvpMJviro/xabwm5O4v9U6A1DDAroWfjVVx1sDmt9QNk4ajUajSVxMm18Tq8ObqTixk9rAebsdt8em+TUEQ0+687nNbWlZwcfUJbwlkY+1KRkBs4eMnND3u7axpprWml9NneCTG+Glk9VyPIRfd4qaJjMzAoE2gdBoNJpEZssceLIrfHO/8tewanut48T2n2G8g4lDZeHQzH3ZrQ3Nr6E19Pth8VtqW0q6swbRW6SiRVj3ffcXmPYwvH0B/LsTFB6M+GfWOEvfhZK8yMPFleZDieGQntYQ8rZB6VHnsvMnqWex8uPY1DUaFrwUWJ4wCF4coRz7agAt/GqqxvL3A8ueOAi/AOc+H7xeuC8+19FoNBpN9VnzJRTuDwg1g68K7Os2GkY9CONs4aWb94Tht6tQV0OuUfuv+briuVv2gbEvwnkvBNv8mjOQzbpVYvZgxMi1a5ZnPwWbpkHRYTiys+JxicKe5erbej8rw9SyH3MaNO+hlgv3O5f95h71/e1fql6/qtK0q/o+7iYVznTvytBpmmOMFn411Sceml+AAZcEryfq1I1Go9FoAjac0q9CVjY9JrAvozGccg/0uxBaDwhsv/htOP1hOO3vKlxXvwuh0wjn8w+8XE3ju1Iswq8xLnT7rSH8Opk9FCub4C4jQ9c9kZNfeIuhQXNo0jmy8secqr6PuxFa9A6cozJq4/d7S2HA5XDmE4FIFjVk/qCFX031seYajydOKS81Go1GkxhYUwinRDhFX5WZQ6vwa44LnnQ1FjmFNYskUkIiZ0SLNNKDE6apRDihsjZsoL1FAfnB/K4hxzct/GqqT6SdXHXRml+NRqNJXKzaw0iVIlUR6lwp4DM1v8WB83gy1LS53x9cPpJICYmu+a2qkilSobI2lEtWJ76UCIX0GKGFX031qSnN79afYOvcmrmWRqPRJBJ+X3DK2niye3lAuIyUgn3BDsoRRyaowvjhSgkkQDIdu1LSAueyan/9ftg8o+4Kv6WFsGFq1TW/5nF7V1ZeTvpU9rvNM5Q5wtZ5sGdF1a5ZVqTqHMrJzlrO1Pybz67kSOjyMUQLv5rqU1lsxurSw5LR5seH4fUzIHd7/K6n0Wg0icjsp1WEnV1L4nudfWvgpZNg+iPRHffVn5XjmEl2y8iO82RGdx0IxH8/vAV2LDTO0yAQKcB6j7bPh6MHAAldTwt9zkQVfmf+Wwn64YRXKz3OUN/NukGW8Rxm/8e5bGp2YPnNc+Ct82DOM/D6GJh4YtXqvPR/MPlClUY6FH6fitVsCudmOzCfZ5zRwq+majTvpZwT/rwG0rLDl68qF70Bd6wMdlQoDpH2UqPRaOorpkCXF+eoBGZUgO2/RHdcgS2awCWTwx9z4etVS2nb6ST1XVoYEHjbDrFst2gczWQPZzyhnOnu2xrscGeSqMJvVbT9Q65Rv7PpMUoANu+LE11OqbjNjC5RVcwxuuhQ6DL2JFlth6jvUKmvY4wWfjVVREKrftCwTXwvk5IGOe0DIVFA2/5qNBpNvDDDZFVXGLQmlAhFo/ZVO7c5Re73WpItpQWm0K12o+Zyg2bqOyMnWNtpkqjCb1UQInhGtmlX8IUYN73FkN2m4jaTaM1fILJ7abXVBovJSs2M71r41VQNX1lw9p54Y7V3cvLm1Wg0mmTAmuY3Hpj9etTCYBWSE7jc0R8DgTr6vMEaxHKnKYsAVb7fYlvsdA/rk/Brx5MRWqj0llQ0XQy6f1VwQDOPqSyCRrnwmxaoY1WvVwW08KupGn5vcLrKeBMk/NaxtJQajUYTK+KdActliAX+ssrL2YkmOdZX0CcAACAASURBVEF1f4NVQC8rUuvuFGeHNzOKgTUqkdP1E1X4jcXzTkkLHc2hrAjS7cKvg+Y8GszIEpVFkCgPUWc8F3dq8LFxRgu/muiREnK31qLmN4ppkdJClbkmnNepRqPRJAJSwvTH4NCvgW371sI6I+vZhu/je30zTJhpa1qUC5/eDEvecS5/9BB89Pv4O+JZMceeb++D9d8F1s1xYsEkWPWpStn7i5FRzqr5dbIz3rNCOYUt/zB+9bazexl8eK26v5/dosYrO9vnq+/UCMxIQpGSrl5mvrgj2Izh8FbYuRDSbamm968LLH91F8x9TqVB/ug6dYxJ/h54+3xVf6uQvtRoKxumwvtXqs+LI2CPxWnPrpEXQj3HWU8El4sTWvjVRI/pEFFDhumAyvhT/mYYhdnDvP/Czy/A/BfjUy+NRqOJJblbYebj8N7lgW0bvgss5+2I7/XtGtBdi2HZu/D5rc7lpz2khEwrV37sXNbkvOeh62iVsrgqmMLuriWwb1VAO5nZRGV6O7hJCb2/vAqHNkP3McEC3in3VjxnWZH6LZ9cX7U6VYXVU2DVJ+r+Lp0Me1dVLJPdSn2fO6Hq1+l0IuR0gEWvqwgZJpt+VN/thkBni+Nbo3aWOn4G3z8I3/0FVn6kQqGZrPxYnWPZu8Fh0UyHt7Qs9Sw2/qiiVbw8KlDGnCkwx3WA/kZWV2vIvDihhV9N9JjCZ6cqhkGpCh2Hw62G93E0ml9z+sse9Fyj0WgSEdNO0qoFNDOntR4Y/+l5+/nLp6FDTL9b7Tpb9oXxeUqwrYzWA5SAXNUY8fZZRzOdr8sNV3wA7Y9XAnFZkRo7Ln8/2L6443D4ky2KQq1kOLNd0+nZ+v0w4DLoO67q1+l0Ipz+aMVrmstDr4Orp6hnNz4PbpkHty221aMs+BuCx2JpjLGmBvjke9V5bpkHHU9Q26ymMabyTFiey9j/KnPKGngWWvjVRI89RElNUe7MEIXm1+yYXbqpazSaOoB0eFE3+zFPRvzT8NoFsHCCiFVzV1MJj9w24dcuDKekq3p7S0LXyT5+1YZpXCTCb3Wyu1lJqSQShtNYHmp8Dxf9wWy/4ULYmQopu9NjSnqNRHwIKxEIIdoLIaYLIdYIIVYJIf5kbB8vhNgphFhqfM4Mdy5NPaE8tExNC79VCIXi9Hap0dRjdJ9dx3Hq3/xeQChBM+6aX5twHa6/tQo5VkE4ntiFXbvztccQoLzFAYcqO9bxKzWrduLHRyT8lgQ761UVxzBwlSiyQt03q/bWGjXDbDfmbwgXycMUkoVNDPWk10iq5Ug8lrzAXVLKxUKIbGCREGKqse8ZKeVT8ateAuL3A7LqIVrqA7Wl+TX/jMVH1Nun/e3fiXLNbxI/r/qE32dEGknRzzQ0us+ua/jKAkKkdeD3lijhwFditPkU8MZZSLNHebDPtEmpPuZsWk06PptUEH4dtIdlRaq/iETzm9qgYkIGn1edN5ah5fw+Q+gTavyyC3lOWn1vUWw1v9ZoCmVF6oXFaWY01DX9ZYHfYdUCl5sYmsKv5Rk5RawoV0zZrl1Dmt+wrVZKuRvYbSznCyHWAG3jXbFaY+3X8N5lcNe6gKG5lecGKe/Wu9eHfjOq75hvjp4aFn5NrcLMx2HeC3D7YshqUfkxP/9XfX//IPQ6Fxp3jG8dNfHh4RbQ/2Ll6V6wFzKbwR3L1aClCSLp+uy6zrL34NObVAay42+Cr+9S23O3wiO2/i1vOxxYDz+Mh9Hj41OfrXMDy88NVRnCTMY3gn4XK+emP/4MgN+VWj6FXNigAzv25FNc5qPM56fU56fMJynz+oPXfca6N3jd65f4/RK/lPj84Jfmsvr2+8EnJU2Kt/MXS5V/2Z7Py28tLC977YFDnFy4HYAf1x/k1Vd+LpfZASQS6Ze8b56gYK/6GEx7/EJOLp7O9PTT+G/2n5S1s5RIDNnfsH82z6m2BwQ8s4xV5suUR/k8/1J1zxBcnfUSrxZ8jlVX/rdPlvKzJ7XculpIP1N9pby5cC9vr5gZ9tHJSsKidfFt4WVgwXuP8kCmEv1uKd7Eb30pnPv0jArlXdLHDw7neXf2Ss6Z1o4sgs1E7ntvAb+keGkgC/gC+O/MLXw0T533gWIfvzELjm/EpZmTmFD0AC2AOz9cwdKUwAvGm4V+Nq3ayj83BdfphpO6cPnxHSr7+VER1SubEKITMAiYD4wAbhVCXAUsRGkaDjsccyNwI0CHDrGreNxYaIRF2bOiovArZcBTsuhw8gq/ZZXYCcUT6xt4aT4c2Rle+E3NVmUBti/Qwm9dxVcCS95Wyw3bwZEd6iVUC7+VkhR9dl3HDCt1YL36tmrmXJ5gTeyRXep72XsxEX6llOSXeDmQX8KBglLyispot7+YXmaBgxtYUtaBQdaDVnwAwMB/fk9RqY8/sIU7DEli5JJT2L9kVpXqIgSkuAQuoT5ul8AlML4FLpfALdQ2l8ggK/U6bi9V4/WREj/bDh0tP26XbFp+3r00pbjMjzCuIRDlF3w650HSKOHW3CeD6jKobAkevJxe/B3vtLhbFS8/HoQoP4sxLAnLPmO7uc1Yb16aB8ZQ5EIypFE+RUXZpPoCf8F2OR66ZWWVH+/xl0A+NMrKpkcTh6x0jjcyxGbZFzZAgzQPPVs3BKDFXvAXpNGrdUPHwz48fCv9CufS82jA+a1LWj5ZDvbRnZumcTSrEQ28EjZB85ws+jRWETZml/2B31iiRIxptJ0WRw8C0Kl5Q3yZgUgcri3pNPNAv7bB4deaZ8fWnjxi4VcIkQV8DNwhpTwihHgReBj10vMw8DTwe/txUspJwCSAoUOHxjk6dwxxeoOyquJrwCYlYbFnZqktIgmGbQ3HFm3Qdk1i0ry7En51spNKSbo+u65itmOzP7O264zGULgvsO6LvA/z+yX78kvYcfgoO3OL2HFYffbnF7O/oNQQeEso8QY72N2XspteFskgvyCfza5OdPFvCSp3bv/WZKSmcOzelmCEfr33wpPJTE0h3eMiNcWFx60+qW4XnhQRWHa78LgFnpTAutsVpXmB72R4WAm/p/Vuy2nnnxzYt2IPfPw/AC47czSX9Rke4iRGFILxFuG3QQualBaAMXS89fvjoqtXKPamgyXi5h2ndIBPZPl1AG4a0QH6DglsOHoInoCxxx7D2GGDq1+H14bTx+XmhcuNc32aBVuyA+sVGKxi+T7dw1gXHN/GA5YwwBx/M8yfyM0ndoRegyB/LzwNFx3biYuOtbw2vXESbJkNwDXDO6meCfjTb3pCB0u5V5rQLi2NCZcFvXLFnIiEXyGEB1XVyVLKTwCklHst+18G4h+Yrbbx2myxkpVy4beWNd+RCD9WG6pEzeCjiQ4zG5EWfkOi++w6hNcm9FrHloycYOHXoQ/LKypj474CNu0rYOP+AjbuK2Dz/gJ25hZR5gt+d2naIJUWDdNplpXKMc0a0Cw7jeZZaTTLTqVpgzQaZ6bSYdFMsOSrOLlzNhT7YOeWoHP98+zuSgHyU+ty4feioe2reheix+rE7GTz67QcCfZ7Hivs/ZW32MHhLYSzYawUTSlpUFoQXIdw57buT82Ckvzg/ebsm9k2nWx+7eex4mTzWwNZ3sIKv0IIAbwKrJFS/seyvbVhWwZwPhD/lBw1iZORe3XzXdcXEkXzG8kLiDVskBZ+6yb2WRgzD30yv4BWQtL22XUVsx3bhWComNXLmMkqKPFyzzuLWL4jj525AaVMaoqLLs0a0KdNI87o15q2ORm0a5xBu8aZtM3JICM1AifRFTZhJFS0AVNwqsk091ZcLiU4SX/FsFrVEX7tqX5jhb2/chR+Q4SZi5WJoScDjh4IrlM43x3rs3d7oORI8H6zjZYLv8bshFP4OScqCL9pSuMdZyLR/I4AfgesEEIsNbb9BbhMCDEQNYW2BbgpLjWsCYoOw9qvABEId7JhqnoLy9+tnAyadAHDRgVQ5Ru1h6zmtVLlWsFXpjKv5G5T6zVt82tn83QVPNuauWfzDGjaDRoZ/j1Ws4d138LQCrO8mnhxYKPyoG4fZtpw/zooKVD/r8J90OmkYNts+1Sv+bwtDiqO7F0FBzZAh2Hq3J5Mtb1J5+h+R92j/vfZdYWDm1RGzA7D1JhSWqCcNdsfF1AemFm91n6p+q/i3PLDy1xpOImWWaX72bhrP4M6tODKYR3p1iKLri2yaN8kM3rzATv+MiW4mMLM1jlwzGkVy22eqX6PPUpCTWIqN+yClqc6wq8t1W+kkYXCkb8neP0bhyxzR3YGlosOq+xqEDvn8pQ05c/062zVJveuhKyW4Y8xcac6aH4N4XfzdOh3oer3oeILiXV91aeBZbtyIyVdZRX0++Manz+SaA9zcDah/jr21akl5jwDP/1f8LYFL6lPKGY/peLdnf5wfOuWSPw6Cz68JrBe09EeANxpkN5Q5ZufP1G9gAw30m6WFsJb50HbIXDDj4bLrUXzu+E7KDwIDZo6n1sTW14/Qwmz48OEZnr51OCpuD7nw0VvBNbtttqtB6jv3cuh51mhz/vWeYFU3FbC1aeOkxR9dl3hzXOUQHPjDJh8YWD7kGvhnGfVsunEBqrNWli1bS8DjfF/h6cjh5ofT/9dyuFs6ph86D8m9nX2eyEtWwlfJk6OpR/8zthnOGL1qOS/GG/szunZrUPvc8J0ooVgzSjA9p9jk83Umv4Xgu9vWiMoyYODGwPbFr2honoANIxRsBZT0HzzbDjpLtX23GFmcIVQzzijsVq3C7+mAnDJO3DWMzDXSMNsF6r3WCaarOmL7UldzPU9y6BN/Ox+ayFAXwJSeDB8md99Ck27qobiSoEXjq3YCOo79t9bG5rfv+5Gaehz4YnOwVMwpgfqzkXq2xqM26TkiBZ+a4pI7easgi9UbGf259jLEBDCxd90Enw1mprE1OQV2Nri3oAgIJEc6nIeqdvmkO09yAZ/WzbKtpzhXsD+3tey6PjL6N8M2qVn086dCgfvhheOqzj9HCv8XjVL0rJvuYMSYx5XL6T5e5QQ9+4lgfJmNJ0LKlEWxZsT/xy83vQYuGez6iMym4Q//o4VyqdH+mHmv2GXxeg51uP83w7AvjXw0klq/e+GEPxUt2ATANOp/u6NsZthHvaHgDb5wAY1u90jgheou9aquk0cAYW2+2EqIyBwD5v3hC6nBJfzGWYfJ96pFI7l2DS/w/4A676Ou3ylhV8gZM5yK407QY4l7I+nQfLZ/Vp/r3CHT18YD0zHhswmgRSWJvZg7E4RObSdaM3j90WXkMLu7GBPp+lyqem3ZI64oqlb+Gz9jpRsO3iUjxbv4IbCQqasL6a/uxlDxEGyWnfltMaNYS38pk8b6GwTfEyNWrzGH/P/GuTolKm2NWobWiipjWQXJk4vwtEoOVyugHbbboMaq37GW6w0qG5PsFbUnNrPyAl+pj7D/CSWppXW8KxCADKgua+MNMO0wZ1asd15LLMCZcXqfjVsE/pczXoEr9s1v6Z5Wpyd3rTwGyl2g39PenILv7Vt7wuqc7b+QSo4FDgIunYBWRN/vMXRxeO1/6+cNPg1lAVIo4kJtra6aX8Bo5+aDsCtaWWc0LMdXb1FsGUdrZs0qtzW0RRg4ib8mhkULeKBtb8PpfSorynkY9XPWCMrODmL26Mc+L2xdya0PkczskQ0jutO9bG2B29x+HTMaXYnTgebX/NccSR+1sR1CaeUgnbsDSQlLfkGX+vvre1ID1BR82t/Q3f68yTbM0sEornnaQ0rlg8p/OoXGU3dwFca3FYLS/3cNqorc+8bRSpl9GzbnJRUQ+MVLnmSOxUQ8evLwgm/oWZx6mu68VgJYd6SwLN1esb28czvi702PUj49YauSyicXnysz91bEj58WgXBOJTwG9+xWgu/G6bC8vfCl6sg/KarVKtf36s8ehMRKeGLP8GaL6p/rl9nBXunJorwu/5bWPi6emN+7/Lg/VbnAZOv74H/9IGXTlZ2V6Bsn/53KXx1N/zvEsjdHv+6JxPL3lXf239R6VGf6Wez+bKQngN7VwRrA5xC1Jma3zVfwn+Hw/THgvc7ZCDSaKrFrqXwzjjYYCR93ThNrVvtQ60UBGzed375r6Bd/VnPn3ffQ+vPLlYbPOnBWsHK7NmFUALHvP/C939zTsgUDVKqvm/KbbD+O+WJL9zBQo21PqEEMru5QH3Bmu65qnz/ICydrCJHgfECYyMlHTZNU5lIwXgJifELhVX43fGLsS2KsTycqeOXd8DBDZUL1PZZDXv7NR3pl79PPNFmD6s+C15vPwyyW8Lqz9X66IdUuDPTDsXEV6Ya54KXYPt8uCl83u0ap+SI8hhd/Bb8o0IW0+h485zg9UG/q975YsGgK2He87D0fyp0UJ5FaPX7AoNPh+HqmZYUqPBYR3aoz4bvoUUvFV5o/TeBY/tfDDk1GKy9vvPrbDjhj6otAuRtU8sn3lmxbPPuan9xXiCer6n5zWgMvzWEiJQ0pWH48Gr1P9y3CkY9EDiPmSpWo4kVG6fCxh9UqLJuo2HFR2q9wwmOXunLfpmB6QqU784B767gAvl7VJvueCJ0PkX5lBTsU1ETWg9QQmaoaCZ+r3I0mztBee1nVCM2rd8Lv7yslguNSAcDLlG2nOY4aMUq/DZoDk2OgeY9wjugxoNT7o+9ImbEHSr602l/h2n/rGiTWhXmPhe8LgQMvx16nBnY1u9CFVZu/XdqPDNDzsUS670yQ7d2OCHy4/tdBAjlA9WyNxQfUbN1vc5V/4WtP6lyTm3hsvdVO2vSJbCtVT8Vcs1KlhGdwx51I8Zo4ddbpP68twdyV7NzceBP3+5Y6DSi4nHHjILdRghNGYHZRG1gThvE4s9r59S/xv6c0XLKvepZHdnpHCzc1BiOexkatVPLBzfBc0YqR9NGzT6trk0jYoMZK9R8Ntb/iWnbZjc56nGm6kStz8CM83vef6GnMVh4DPs4q1bYGo/T3h4G/Q6WvF2936NJbsr7U6PNmiH4bG14474CHv5yNZ6Na3klFeaN/oRhI05VAsF4SwzZM/6txhGTtoOh77jA+gWTQtel62j1P3G4ftRY/0PFeSpO+ol3quVv7qlY3mr3eeKd6sW2trC+8MaKzCaBcIgrPq6+2YPdYdfEHiZ16LXwzX3BmdJi7VRu99VJyVAvLpFy3A3qY+eSt4PH1g4O6aRb9oazn4EjuwPbxr3mEA84BXqeDYd+jbxeVUALv1Y7HBPreqg3r9pO7RsJyeANb2oA7Z6hVuE3lO2aOYjZ71My3Ld447Pcf0fba4dsVmCx97I8A1P4tXaSdvs48xh3duhrWusWi6D1muTC7BfMtmW+NBvtvMzn56WZm5gwbSNpHhfPDGkNK+CEri2dNWHR2FrasfZj1bV9twq/RbmBaedQY5x1Kt5p+r4+EQvH9miOtyYXiYvNr4PvUqywtufKcgAEjcchrl8DPh16BCgrqrxBhBR+E8DmNRzJoMH0ZIROExlO+C0XziKIEqGJDmtop8qEX/tLi8fB2cEUMioIvw7PLS3b+bxB5SxCskYTKeVpiI3v8tkLL6t25XHPh8tZvfsIZ/Vrzfhz+9B861ewgviMIUHCbzX7qyDNb24gTFUoraP199RmeLOaIBZRZaI53pUS0OTHw+bX/hIWy6hNkaaUtv6mUC+ANRDNp5633AgoyXeYCrCGdQnVcSVAqK9wVPAcrYfeuClp6hkW2uyDSgqUFgOCf7d1wPH7lFYx32aL5y1SL0WeDJU1rqxYGem7UgLClaZyrMKnRUgox1usHB2KbLbo5v/KmviiXPhNDS5ntwkzNXNSVp5gwyokazTh8Puh7GjwbEVZUXm82xXbD3LB9J/IyUxl4pVDGNPXsFk023BI4bcaY4j12OrOVFnNJo7sDNhkhrLhtf6e+q75NceX6hCV5tcd6CfLjsY/fFwss7Rax9bKXuwi0vymKZt4bymkGCmVhSu6kJlhqKfumVGwYwEVMoFab3CoqR/r4BnrWHyxYtvPgeVPbqz6eY7WYu72sAglQH1q+33P9Ibpj6hl65/N+qbp98KkUSoto5WZT8CjrWD/eniyGzzZBf7dCR5rVzENtsaZoOQjppbXNkiv/BieHxK8zfxfvXwqfHO/Wl5g2D5an11KWsXzPdsXts5TEU6m3BbYLlzBThbarEUTDVNuhcfaBqL6FOepvmCbigIwf9M+Tu7WnO/vODkg+B7eGmiDoZQOqVnO2yPBOv5sma1sic0oAdFij6Zy1Jbx1O7sbdUI14aTW02SklH9ZAtLJ0detqxIOYXNe0FFooqHv46VeGl+K0ucESQkVyJfSR+8cSYc2gyPd4TZ/4lNPQ2SW/g1Q2zYs5GkN4LL3oNxr0Kzbs7H9j4PzjFyWDfvGb86VgtLCJHDW6p+GkvIHrJawq0Lq36uWNO0a/B6i94VywRN07nhyo/Vsq9MhdUyOetp9W0Ka3tXQFlh8LmqOsAkG+Y9dKdVnC7+/Xfqe+8q9W19hlbP4/kvqm8zhFIrSxpNT4YSQuxs+B5WfhJY//33cOMMGPEn6HOB2uYUOk2jCYUZqu+w4YBTsC+oDZ3QqSGvXD2Uxg0sWtAjltkka//zp+WB5epElBlxuxqDQIX7A9j0Y9XOZdrUm1iF9Wu+huumVtxvZjut7yZipk9JdTD7qeumwh0rKy9r2rkueBkymqhoGrHG7H8htsKvtd20Gxq6XEoaXP4hXPGR0uo6MewWNatweCvk7VSCcHUimjhVN6Znq2uYf9wWDsJrjzNU6JFQb7ZpWTDkahUpwikIfyJgdtAdT6xeJ2U1PG9yTOgXgtqgccfg9X4XVixj18x3Ha3eOK1CUKeT4Njr1YuPiWk2YaW+azpihTlgpDcKtB9vsfICbtVfrRcb9/fEPweOc3ugi8UD3jyu7ZDg+JApaYHnc8ypznVwpUCH442wUe5AWKHqesdrkhOjvfmOBpvq9GnVAGHvF6wCk1VTau+vqkrDNjDMiLIQK4e37mPUd+POgX2dRkCrvhWP6TpafSdqpKNY4eRYGy3eEhXSrv1x0b3w+EpUKLBY02FYoI+Nl/lmOBPL7qdDt9+E3p/dUoUyRQZkl/bHx6x6kPTCr9GoqxO5IRZ/jnhhdmppWdWrYyK/3dufndNUotMf0epYAAGh1toZOGkWNZFhtreMnGDNb0pa4B6b99dud2bXzHpLnEP0mIN+ukUjUNnLidkOtOZXExVGmzLam9vn4Fxrx9pnxssprNw5tJrKF7MfjKaepi2qP87T8rVNrKI9RCtjCOEciSpWmOdNZMd9e6jMGAvqWviF6jWAWPw54oU5nZXaoHp1DLKRrGY2oVhjf3ZOwq+TQORyB+J0QsAEJkj4ddD8aiLDtJNLz1HtR0olPHgylAbXnRrQ3No7NbtmNlxEFut0mJSEbKPm4O4vc96v0TgSps9zmkmwamPj5WjsFBawKjhFxQmH+ZuSQvNb3WgPYdL9OiGlc78XK8zzxku4jgWuFBWWUgu/ccD0NK+u1+3+dbB6SmzqFEvMTjk1S2U/Ky2EX15RqWWjcWLbPD0+9YsF9meXmulczo4rxTktqfV8Ts5t1U0lmgz4faqdgSGYSpViO3+vJYVrukorDRU7eKsmbcdCQ/Nr66Stz8mq+d23xhIpwvbSUy78as2vJnLC/uP3rlKZQtd8oQbrxW/D8g8D++Ol+TX/N9VNx14V4dfU/MbbIau2cXKs3bdWpZY2fRYqY/13qk+KVsY4/Kt6sYiXWYJ53oTW/LpVFsNZT6n1GNc1uYXf3cvUd3XCZ6SkKcHyg98524jWJn4vIAKd28LX4Ku74IfxMPPfkZ9n2XuB5f4Xx7KG1adR28ByZlNo1CE4PEyjEDZWjTsFcpsDDLlGfefvdiodIFG1/InE7qWwynA6MzPrffA7FX6syTFqvUkXlcbY5VHPrO1Q6GhkUrQKp1/e4aw5aWqcR7iD02Na01Tb0yebtpfa5lcTBcWuwGyStIf2Em7V3j+8Gt6/UkVemHIrrPsqUMb+4tayHzTrXv2KZTZTNvVmTO2qCgfmTEj749T3gMvCH9PLSHff+eSqXbOukJKh7o+1z5j2T/juAZj698qPlRLeuxz2rgz0V+HoMtK2IU4+JmZ9mkRYr0hp0CJ2bcL01dm/Rn1nNo3NeQ2SO86vaRbQdnDVz2Ht2MqOxtwjsVr4vept/phTVbiV/D2BfdGkDpQS+l8K5z4X2juztmjYBv66R0UEEC4l4DywHRDqzVGEeL/7/XdKQ+hyB4cNOu4GmP10YL3fxSo9ss8Lr5+R2PbPiYJp8nDZ++r/sPA1tT78djj1QbV8/TR1/92pSlt/3feB403h152qZit8ZRWn5wZergZgM/bjvb/CExZHnSs/DjjlmGibX02UTJy5iQu8LlY3HsOgGyfiysiBh4w+/qop0OlEmHxhINKCNUxYp5Pgyk8qxor/w5zYVC4tC+7eoJQuT3evuhbW/L826xZI6xuOjidEXrYuY75QeIsDSrJSI+5vSYHzMSZeI/36KfepTyRc+YkSrn96Vq13GhF9nSPh5HvguBshrWFsz3vPhtidyzoTMeIO1d5jSHILv7GwJbG+bSeaVtAUfk3Boao2rH6vEhATTfA1sQtGkWjy3SnOLyr2oOJmGXeK+vMVH6laHZMJq6Ol9f+RmhWwv7bff6tdpKlladBCvWz4Sp21WtaXFrutd1ojKqDNHjRR8L/523j8m7VclellUPeOuDIbBxfIbKLarbXtWZOzZDSOf5+ZkgZZLdRyVV/MY+H4XV8pt6suCYwr5stCuPHea/F7iNTu2+WGBs0qXj8epDv0kYmEVfiNQ13Dmj0IIdoLIaYLIdYIIVYJIf5kbG8ihJgqhNhgfDcOd66Ew+wsYpVpJ9G0gmZucFNwsEYviCZkl6+s/mfysWP+8VxhUupqKlJuQ+gJHlBDZUsMdXxGjpG62iHagx27saJ1dAAAIABJREFUvaKTsJwkwm+97rNriO9W7eGvn63g1J4tyHB5cTllwjLbtvXlu9QSF7ymsoAKUb2oQ7Fw/K6vlEfUcEjaE24sMPdHm0Ut0jTB9Z2gzKyxvw+R2Px6gbuklL2AYcAfhRC9gfuBaVLKbsA0Y71uYXrJVufGWht2omWO8pepBmR20lab5Ggct/xlkQsu9QWnNpGSXn3P6mTAGjrJOqBG6lBj2iCmN1JalrKiCIRfV7CJi5MXs3l9X/0WfqnPfXYNsHJnHne8t5QB7XL47+WDEKFCVZU7b1rauHUqPJapY8ORklb1TGRx8qavF5j3xDq2lwu/YcaCsirKF9byNdmGEo2gzKyxvw9hRyMp5W5gt7GcL4RYA7QFzgNGGsXeBGYAERq2JAibZ6pvdzXSE9ujA1z8ZvXqFEt+eRXSGwbquHNRYN+G7yqWP3oIPr4eRj4A7Y8NbPeVBml+fX5Jmc9vfCRen59Snx+vT5ZvK/P58fr9+PyqvJQSvwSflPilxO9X69Zln1TlfCH2IWW557VaVWvWbea6tAn3UoJ5dHA5c1ktDNmax/FAkd9NBrBkRx5zpm1AAqP3l9L74EY++WgyF6y8hdz0dnzR8wkONEicpB+xiEUhpJdz1t5PVul+jqS1Yn6733Pq5idZ1moca1ucEfb4Lge3MhZ455edFHqKucnYPmfdLhYUrgt7/OWl6bQCNh5x09XIsDd3awE/f1/5sbfjJgVl9/jy3J3kpwfbQLbM38kVAJPHsfeKabTsVkkWojpMve6zoyV/L3x5J5z/YkRTp/uOFHPDWwtpnOlh0lVDSN9rRISxCrjpjdQsWrnHvGUMmPGvwHJNCpPuVFjwknK+O/2RgPNaZaz5Er65F47sVOvJLGiFwnzuzw2GbqerDJImh7eotNLHXh/IDmriK4MJA4PPEfE1teYXCBZ+43AfolLnCSE6AYOA+UBLo5NFSrlbCNEixDE3AjcCdOjQoTp1jT3mza1O1q6WfQLLB9ZXrz6xRvqQadkczOxMerex+Ipyydy3BE+pMn948ce15JZIjpb4KCrz0fHwz9y2axqLdhTwcMN/UFTqo6S0lBnSz/MztzLhx2/w+vz463G0rwYcx/0pm3i/ZCS3pEzhuU19WLNRPdcVru68nArHrfg7CMgp3sHaX35gcj0LHtCSQ9yZPpMjMoNWBav5Zm8j2qYsY0uulwkruoY9/reuHYxNhXcW7GSdFFyZlkYDUcK2X9fz3IaNYY//SVzMZe4fmbFvAM97ZgHw0qbGzNpY+bE3p7pIMf7Kj8/Lx09waup0JFkpw2lIIc0LfbQMW5O6T73rs6Nl9lMq8sKy9+D4myotWlzm44a3FpJXVMZHNw+nRXY6rDJSuVuFyXGvqXOa3ue9x8KCSYH9GY1V9BIz/XBNkJYNhfth+3z4dWaEwu+UgOALyS1ohaK9JZKMVfAdfhvMfU4t//JKReF3ryWNsTVleyR0PEE586bnQFar6I6tT3Q9DeZOUKZEHWPv+Bex8CuEyAI+Bu6QUh6pkM4xBFLKScAkgKFDhyaW2OQtVh651cF6fA07vBWUeNl+6Kj6HC5iX34x+/NL2J9fwsEjR/ka+M+B43nuiXlAIETZDe4v+avnfzz3/Sq87kyy0lPI8LhJEUpY8MhSGmZ4aNkwjWx3KqyHnu2acG27TqS6XXjcLlLcglS3ixSXwJOitnncwvhWyykutV8IgUuA27bsEsbHRWBZqGVVNlBOGNshEPxFCGFZNr4pX7Bsq7y8vSkLcT6XIIA/8hvLOQRj4OH/0LYBcFRtf+Tcnjx6/FlVfYSJSe42eBYadhwI2+bxp+HNYQGM6N6SLb+L4LeuLIaP4Ns7R6nU4Yufgim3cfmQNlw+NpJ7dRZwJxcWHoAn1QDz5t//GF5z91galJRAn/PZdNE5IQpdEMH16wf1ss+OlnJTl/CJTf719RqW7chj4pVD6N3G8II3p67bDgkU7DZafUw6jYCL3oAPr1Hhx+7dFJOqR8WAy2D6o2o5UvMHe8g/LfxWpGFr5+2nPwJHD8PSd9S6lMFKNOvMY8M20V2zUTu45J3ojqmPtOgFd8dPoRiR8CuE8KA60clSSiOAJ3uFEK0NDUJrYF+8Khk3vMWqs6oOcXZ4k1KyM7eI9XvzWbengPV789m8v4Dth4s4VBic1jLV7aJ5dhrNstPonOOGPBjatTUPde9DowwPDTNSaJThof2GTTAHFt9/Iuk5Fv3XmqPwPvRv14i3rjI0B8VH4HEY3ac9o4f3ivnvq3O4U5UNoIGojzFjTYcw07vZtBWP2GbXli61qo5m1unCSAZm00HCHrEjCam3fXa0RBje7qvlu3lr3lZuOKkzY/patG1mn+4OM3Vtb+s1jdU5KFIljL2cFn6jw3rPfWWJGw1J40jYf6pQ6oJXgTVSyv9Ydk0BrgYeN74/j0sN44m3pPp2TtbjY6D5PVBQwpJtuSzZdpjF2w6zcucRCkoCHXfrRukc0zyL3/ZpRPsmGXRokkn7xpm0b5JJ40wP5dqdwgPwJJzSuwOnHNcp+CIHVYipdGHThjjVvzzmajXsousTKRlQZplOr4+pcq2ZASE4SkhEx5ttxiYQRJsKNShSRAQDixmZI17pZOsI9brPjhazTVTyP91yoJD7Pl7OoA453DumZ/BOb7ESfF1hfMPN69RWP2kVuiNVwtjLaeE3OqyaXm9xsPBbH5Ui9YxIXlNHAL8DVgghlhrb/oLqQD8QQlwHbAMuik8V40RZERzcFGyzWxWsHcbRg1CwTw3UESa7KCzx8vPmg8zecIBZG/azeb8SrFJcgt5tGnL+oLb0bJ1Nj5bZdGuZTaOMCDvXI7uM+jloLJw8WCE4CUb5eQybMC38KlLSoMQiDEaqzfR5oehQICZnZRQeUIOZmQCiJvD71G9JSQtMEZvC7z4jw46vTLVvV4qKcRoKsx3ZtWERTD0HYY0wEsmUvRntIVRik+ShfvbZVaF81sFZGCnz+bn9vSWkCy/PXdIPj9toO1KqmL3710XW99W65tcq/EYYkabIluI+nICvCcZq2lCSrxzDzRi9iRbzX1OBSKI9zCF0jr3TYludGmTyRSotZHUznNg1Uk8Znv9/PxRSA3W4sJTvV+/hy+W7+XnzQcp8knSPi+M7N+WSoe0Z3LExfds0IiO1ihqs0qPwkmGL7PT7zG0LX4Mxj6nlQ5vh+7+qZVOLcWQXTDwxeFuyU2ibKY70Df+zP8CKD+C2xZWnuiw8AE8a+5v1gFsXVK2e0WJmqRqfFxDoTQE3b5v6PvxroH0/uM/5xWr99zDtIbXsMQR3M8V0izibzZj1TnKzh3rbZ1eFMCY3L87YxPIdeWxJvxymjIBrv1Y7pv8LZj2hljMiCIdsvqQ271l5uXhRFc2vNfqPJnoaNA8sT74Q9q2Gu9ZBdisdD74OkGTBWy0c2aU0oCOrGepSCLj5J+WM82zfwPayoqB0fGU+Pz+s3st7v2znp40H8PolHZpkcu2IzpzSvTlDOjYm3ROjQducpu4yCno4hKY65lT1bdXE7beEkTIzZxVYBL1kS3IRjr7jYOXHkWt+N05V3/l7wgu/JgfChwWLGWZ6Vgj8pg7DoHkP+PyPaj3P4hleWugs/O5fG1g2hef2x8Lvvw92GoqUG34k4vz2Oe3h6AGtwdIEqMTmd9WuPCZM28C5A9rAOmDrT4Gdy94NLJ/zf+Gv0/54uPR/0O7Y8GXjQZDwG4HW0dRath6gQls27hSXatULblsMyz9Q96rkiEppDTDidji4AVZ/rgRfUBE3slsFtO9X1X/LorpK8gq/3hIlwEQyDR2OVobQ27If7F1hnL8Y0rLYk1fM2z9v4YOFO9ifX0LrRulcf1IXzurXmr5tGxKpB3ZUmH+8/peEMHtIhew2wdNj1ikc0z7O2okmW5KLcIweD6unRD6Vb2rOw01JJkISDVOb7U6D7mMs2y2/NZTQH8retsPxVatLNAJz816wa0nSa341FkJofku9fu76YBmNG6Tyz/P6wBP24yxtqF0EYcNcbuhZi1FfrMJvJNEefIazdK9znRUkmgBNj4FRD1Tcnt4I+l2khF+T8tTHhuY3O8pID5oaI3klGm9xXNM57tx/mOe/283Hi3bg9fsZ2aMFVxzfgZE9WuB2xUHgtVKetrmS3+epJFWvKfwECb9a8xtESroacCLV/Jr3L9KUmLWFlJb0xO7QTjC+UufttWXzCIH2rm1+NSYhsvq9OGMTa/fk88pVQ8nJdOjb4pxdKuZEq/nVWd1igz3zn3lfTX+autB2kpQkFn5LnFNWxoirJs1iu2jDhUPbcfPJx9ChaQ05LkHgD+iU4tXEngveqoE2hZ8gQSzOAntdo1z4jdDm19SchxuYattRwldmEX5T6pbwa7b3JI/2oHHA8pK65UAhL8zYyNn9WzO6d0vw+yuWj3N2qZgTrc2vWUYLZ9XDrmAy76t+uUh4klP4PfQrlObHVPNbXOajoKAEM2rw5Eb/xXXFh7Ro10Vt2DZf2QUJl4oEEc/sP+bUS2W/LyUN1nwBP01Qb6nW7HR+L+xdDcvfD2yLNkxVfcfs1FZ8GHiezXuELm9qfsNNSW7/JXi98CA0aFr1eoJK8brkbSWoHz2oArcP/T1snAZdRwdP2/06EzZNV8uulNDmLis+glPurV69Yk255lcLvxoD01Rn2btw0p+Re1bw9zmNSXW7+NvZvSFvh8r+ZjL7P9BmULCDb10QYKwvfDsWwPYFqk9vfzw066b+r93HBPxQ1kxR33XhtyUydgXT7KeVGeRsI8Kgvr8JS3IKvwteVt8x8sxdtPUQ93y0nNNzB3O/RwmRrYo2wtyH4OI3VaG3x0LZ0cBB46OMnRoN5W+dlWh+zWnAqX+ruM/vhTnPKIcuk9ryYk40ht0Ca79U4Y9K89Vn+qOQuxXOeyH0ceU2v+GE35+D19dMgaHXVq/Oy96FHx8O3lacp55xp5Ngy+zA9skXBpazjQQo7tSKmt7pjzoLv6aGrVeoDGtxpLkRTaJ595q/tiYxMfu50gL4Ty8EMKt4Mn8/uw8tG6bDkyOVk5LJtIcgqyW0GRzw34iHX0asadJZKVbSstV/e8rtsH+N8kM55//g4+tg4BUw9r+q/HQjyk/T8OnKNZXQqJ2KnlRyRK1vmwvf/RXyd6n0xGaiIE3CkZzGcWWFKkzJgEuqdZriMh///GI1F06cR0mZn+FXPwL/yA2EEjtkSXNpFXwh2MEs1kRi8zv8tsDyzT/B3w+runc6SWkIy45Cy75q2z9yK9dqJhNjHoM/LVcDYneLo0jp0dDHQCBWaLgpSb9XeYzf+6tRPgZmEGbbO/uZwLajRoxP00v52OsrHte4s/r+2371svb3w3CpxQve52DvbAq/lb0IxIsBl6g6Dv19zV9bk5g4/H/6tWrAVSd0VCtWwdekrEgJk6nZ8VVSxJK2Q1R4zfu3qdkcM4bv3hUBwSxve6C83/v/7d15fFvVnffxz/ESO4udfXMSh4RAIASSQIAEKEsT1rJ1YQrDtJQyZWYKU1pKp/CUp9PpM/OaaafrlBaasrUFWpaWwlAIgZSyhn1LQhYC2ZzNDonjLN4kneePc28ky5IlS7Kka33fr5deV7r3yvr56urop3PP4j7z9fPyH2t/UjPOHfNv74YveVfM/OEwv/KWmmAVsdKs+e1s67k9bBo2frSff7rnTd7b1sLlJ9Zz03lHMqTKP5xp1BTkYna5pH/ba2zf0yWX2NceMCg6PFRZOXR2uC+Aiupg1Hrkm39MuszulyKpLU9ztAf/3PTHyI2fiCQT/ntZGVML4bdVbt/nlgPjJq1I9N6XlXXt+Bhqg/IhXfeJbS9cCBrmTGIlSH6/c95UKspTnCeRUPBGuDlYLg1yQxH2JNSe9XegeIxxN7/M9qeD1/EtagH7dOdIqC2rtjhPrtzODQ+8Q1mZ4Y4r5rLgyLGZxdBnyW8anRli///Y+2UVrtNTqF3tlVKJPT6pklo/GUw52kMbVNdGa+1zMfpDqN39vdgrAX4zhrD396uHdn1Osvc+drYrbzi/Lgqd/IrESpD8HleXovOxjbgyMKjncEW1a+aRjLXed6CSs5zyy9c2L/kt77vRpCR7Af10Z8lPBnrJWssdL6zn3/+8ilkTh3LL3x7LpBEZjuLQl0NapdPTNPb/75L8VroEJtTWPSGSrmKPYar302/mknK0B+/cNMYVnrkY99f/sRf7Psc3w4mfjjtp8htX8xvPbwqhGQGlGCQ6R9P5DEZCwT2H47/bbNxoFuFOwPbpUJ8lyS8z2/a4clJXoYpaab47axf3ulYzHLH82/++x7//eRXnHT2O+/9hfvLEd+hEt9y+HH52nLvF68vJDF67wy17TH5jfvV3SX7LYfu7sOV11fymEnsMN74Iy34Bty+MtqcF1xb4jrOjHdne92Z629cEvzzNdY7wbX7VdVLxj3u4HV78Kfz5huzi9JPf2CsB8TVD3Wp+k4zrHFvz27i6+/b1z7qlCn4ptK1vwVv3dF9/93mw6Az3uUwk0ulGuglqk6/4y+33fMot1z8Hi2+C+/4m8X6SncqY5Fe16kWvNL+hKqrpzbi14Yjl6w+8zd0vbeCqU6Zwy2XH9jwV8YU/c8nEjIvdlIjjZ8ER53fdp09rfr2/3dOc9OOOhuOuhFO/EZ2XHro20A9am7d8O+YSN8PPtIXu8ZM3QcNr8NG66D57GrqO4OAntk2rYdvbsOyW6LbNr7jlUZ/s+jqr/je7OP3kt+5YOO4Lbl18m8D6+TBoVPTxKdcn/luxye+Bnd236weTFIuG1wFYNnhh1/VD62Hrm10/lx+/GY7+GxjhDU0Z7uheYxoUx/TQkfvVX8GHXses2NkbJXvVw+CUr7mZY4ttGEjppjSzm0goOj93CuGI5YYH3+FPb2/lG2dP55oz0hgaZuJxrgdovP+sh3av93BfTmYQbncFeU81FwMGwQU/6b5+0onRcV8nqSdwjyYcB5++HfZuhx/GjIYR+8Mm/n3229om+vHj7zvtzMTPyZTflKK61g17tPrxaEc336CRcMpXYcnNcOI/Jh9eLVWzh1A7TE7vsyXSp7zOot/dtYAnqp526w75mPvB/5sXovud9M9uHcCGF13NMMCU0/IYbA5NnOs1X0sw9bq/7szvuml7JXeMcdPeSyCUXs1vJOwKgDRqqKy1fPMP7/LwW1u44azD00t8exKbi6Yz/3qmshlJIva4aPaf9MSfS12S37gk1x+9IVniaMq61q4m27c3QnGjm1Qm6BBjTHqTQ8S2g0yWwKstoRQD7/yMVMU06amo7n65vyLus3FwfYDP4/gyJJ4uy0uJK73k9+DUv6kTux8uWctDbzRw3YLDuPbjh2X/2rFD+/Zlza8/tFUmko0CIcl1S35bE9+H7tNfxko2vFyoLbtxoTvjEtKkvcHTeI3YL9VEw7BlOZKKSK5s37WbkC3jM/NjJuiprO6e1Cbr/BvkNrGpRqoIcmIvkgOll/z6Y/Cl+IL+3aubuOWZdVx6/CS+ujAHiW+8xlWuQ8b25dExV5s3u3U7388u2Wlrzrxw6/JFoAIyLfHHqTlmMPmtb3fd1tLg3uP1z0XX+e9105rENTb+0EuZ2r0+7kdNVfdmD+mKjW/3BtfpLfZc7csh/ETSFQ6xa82LhEw5l50cM+Nf/Kgn/rqE9wNc/qWaXEE/UKXElV7y60/nWz0s6S4vrdvJzX9awWmHj+b/XTwTk6tev9M+Hr2/+Juw6HS47RR4425o3wv/M8etu2Wu65CRid0b3TKSYWeN2E5yujSWHmO6nk9P3uSWW95MPH30otPhzV9HH7/3iOsYt+6p5D96Mr1S8NEHsG9H9AcWuPc4tj3ggBq3rK1zy56mPI2drvONu+AXJ8Lb98bEqfGhpfC2PPkjZrS9TahqBEMGD4mWZdXDupf9scP8xY56EvuZCZoJc7s+rp3Y9bF+oEqJK70Ob/40tPE96j3b97Txz797iymjBvPzy4+lMtVMQL1x8W1wxs1uSs22Zpfo/P4y97i12SUk0850SdD+BD3p0+E/b9LxmT1/6unR+0Gu+ci3Ly52NaFP/Avs2eLW7dsR3X7Z/fCXf3fTjR7/9/Da7dFt+5ui06z6HW8Abnjf9c5+7vuZjw6yz5tqc9Zl0XUX3+aGs/P/5rij3fLIC+HKxa7DTDLVQ+Hqv7qmFHd5vcW3r4huz6bJjUiOrFy7jgmAueJP7sfp3z/tfmDWn+h+/F31tJvut6oGppwafeLgUTBgiGsWNDHDMrQYfGoR/OoMVyYdfq57vHcb/PwEt12fUSlxpZf8htpcoZZgHNPOcIRr7nuT1s4w9//dsTHTFedIZTWMmuZuvrJKF9PBRGSmS34zndbWryGMn642XWXl0cJfBWT6xhzpblvfgme/537YxA6VNP0cePkX7v7EE7omv6H26Ps/9qjo+iFjYFi9t0+W58PQCdF1Qyd0fewrK4fJ81P/zbo5Pbyean6lsNbv3M/Wnc20VdUwuG6GWzluprv5Jh2fvIKgbg5seD7Yk/wMGuFG69m9AernuZFeqmth6hluqLNUzSJE+rmU1ZrGmDuNMY3GmBUx675jjNlijHnbu53Xt2HmUA8dcn64ZC1vbNzN9z59DNPG1OQnnsqBXvLjJSn+JblMa/r8JCkXnTV0aaz3/HMrNqH1+c1n4o9rqC2mI2bc++Y/zvh8aO8aV1/TaA9Fod+V271w+/MfMrCsk8rqDGff9D+nQZ3kwudPNR7kjnsifSSda/p3A4lGw/6xtXa2d3s8t2H1oSTJ7xsbd7HouQ+49PhJXDCrLn/xVFS5Wt6Dya9X25BxTZ+f7OQgAVENXu8dTH5bu7fT9dvzdpnz3bj9/KHvkvVEz/hKQGvXuPqSP4ygvmyLwd30p3I7TTv3tfPQGw0cPqKC8kx/vPuf02w6HRcDv11/ou+CgP9rItlKeV3fWvucMeaQvg8lTzq790Y/0BHi+gfeoW7YQG4+f0Z+4ymvgjd/Ayv+4B77nS+WP+RmG4ptj5aOzhwkO/4wOUp+e8//onnsemjdlXif2BETyirguf+Ojp+brCf6iz+FuV90X2gbXoCTvuIuY/bEWnj+R13j6gvbl8NLP4O2lr5/LUlLvyu30/SbZRsZHt7JnOYlMGp66if0Z36HPZXjIt1k05vrWmPMu97ltaTz6BpjrjbGvG6Meb2pqSmLl8uRBDW/31+8hk27DvCDS2blvp1vKpFOsGFo9xKH8bNhzAzY+CI8+/3e/71cXOZe+B0YNhmGjM38b5Qqv+Z+5R9hkzdd8fxr3fLUG1yyWzfHTbF6zveitTP+sjauLe7Iae59WPlHeP4H8KdrXLK84fnUsTRvdB3bAGrGZfd/JXLajW7Z0uBmhnvOO1/Hzkz+HCm0lOV20ZXZaWrtCPPbZRv4at17bkVs+/neOPUb7sfo+Fk5i60gjrzAde4bHfMj4GPXuzJowrGFi0ukCGSa/N4KHArMBrYBP0y2o7V2kbV2rrV27ujRozN8uRyK65CzvGEPv162gc/Pm8y8qSPzH8/J10Xv106AEVPgy8tcx4RMLnX7l9qzSX7nXglffReqhmT+N0pV7BfmVUvgO3vg7P9wj6ecCt/+yHVG+dpymPePXacyHjOj67BL4M6HG9bC5JPd+dC6261P59zw9/nMXV2HKMuVM26CQxdEa3wBRh0O0xbk/rUkF9Iqt4uuzE7To+9sYfeBTk6u99r6XnxrZn9o6mnw7Z3ucxpksy6Fb27o2kE1tgwSKWEZJb/W2h3W2rC1NgL8Cjght2H1oVB0KKZwxHLzn5YzcnAVXz+7QJfIks0uVFGd2diuoSRtRyU/kr2fyaTbLjH+fEjn3MjFD6FUKge6Yft8usRatAJdbqdgreU3yzYyfWwNE2u8rzWVgSKSREbJrzFmfMzDTwIrku1bdELtBwvF37+2iXca9nDzJ46ktjrFXOh9pcu88nFzzGeT/KrTUWFUxL2HKfeP3aeH3uWV1V1HfEgr+W1PP45MVVR1HdJNyW/RCnS5ncKbm5pZubWFz82fjPGbtgV9tAYR6TMpG7gaY34HnA6MMsY0AP8KnG6MmY3rM7oB+Ic+jDF3rHWXgisH0tLWyX8/uYZ5U0dw0ew8ju4QL7Z2wsT8Fqmojo4A0Bv+c8q7j2MsedDbmt8u+/TQBbuiumtTh3TOjVx0fkwlfhZADY9XFPpVuZ2G3y7bQE1VBZ+cMwGWtqvWV0R6lM5oD5clWH1HH8TS91q2AhbKB/DLZz+g+UAn//f8GbmbvjgTsTW0NTEVMxXV0LIFfnYcfOFxqBkLe7fDXee6DlBXPtG9ZuOjD7xOR0a1HoXSpelKGl/AsT9Sephym4oq2L0++njJt2DeP7nB6u+9xG3/7D3R7S3b4LcXpx9HpuKTXdX8FoV+VW6nsHNfO48v387fnljP4KoKr1OzrnyJSHKlNcNb80YAWipGcMcL67lodh1H1RV4Fp/JJ8PxX4IVD0U7RgEc+3mXrK99AnZ96JLfj9a5+7s+dAV8fNOGptXR50phVAyAs/4Dwh3pzRAV+x7Gvv/x5nwO3rqn67r2va6D3PtLuu//0ftuOfEEGHt06jgydezn3ZBKeza7H2XHfLbvXkskgftf20xHOMLfzZvsVmiiFRFJobSSX+8y8H0fDiQcsXz9zCIYB3LgMPjED9wtVt1sOOWrLvn1JyqIvdTtNd/owr/MPf+avotXUjvp2vT3HXWYW875XM/DD006Edcm2MKgkXDgI/cl74/lGc8/V875z4RTeefM+FlwwU/67u+L9CAcsdz3yiZOOnQk08Z4o9P0MIuniAhkN85v8HgdgJ5YvZvLT5xM/cgMp7/Ml4Oze3mJTJfe/gmmu83l7G5SXIyJtgk/OAV2W9dzIhLT8UyjfkgJeOmDnWxpbuWyE+qjK0Ptanv+f1I8AAAaPElEQVQuIj0qseTXJQQdpoovn35ogYNJg99uLZQo+U3Q4elgwqP2bv2aPxZwqD35CBAHfwjpXJD+64HXGxg6sJIzZ8RMyNPZqppfEelR/232sGOluzRcVgmDR8GeBg6sXsogYOHRkxhTG4DC0a+1W7vY/Q/bl0e3tWxxEyDE2t/U9XlS/GwPIzx029er2fVrfhtejzabAJf8DvCuZuzZ5JY6F6Sfaj7QwZMrt3PZ8ZOorix3Kw/sgi1vwoQ5PT9ZREpa/0x+9zXCrScTP3SU38jh0pOKoK1vOvwOU+/e726xnvpX+NLSrus2e9PpVhZ5cw6JGnOkW9bPS73v8Mmwe4ObRe2DpfDIl7tuj60Ffvs+t6yuzUmYIsXm0Xe20hGKcMncSdGV914CHXvhwO7CBSYiRa9/Jr+tzYCF+dfCslu6bZ4wqochpYpJ7BSUf/Mb19Fp0Cg3hFXsxAK+imoYMq5vOzhJbtXPg+vehWH1qfe9cjHsaXDTlVYOhBd+1HV7KGYc4PIqGDszvREnRALogdc3M2N8LTMnxJzjW153y84DhQlKRAKhfya/fhKQrPd8ENuDHfKxaDI84TjYtb77PqE2GDap+3opbsMnp7df7Xh3A6if3317fPvfOl36lf5p5dY9rNjSwncumFHoUEQkgPpnhzc/CahKcsk3iMlvl5nDqhJ3eOvUED8lI1Fv9vgObzoXpJ968PUGBpSXcdHsCUn26EVbehEpOf00+fVHPUjy5R/E2c9ie+1XVCcZ6kyDu5eMROd27DjQoVadC9IvdYYjPPrOVhbOGMPwwWriJSK91/+S34Y34MX/cffjJ4EIsrKYt6qiGloa4OXb4PFvwP6dbv3WN1XbVyoSvc/+j74d77mRTsor8xuTSB68sG4nu/Z3cHFsre+2d+DRr0QflyspFpHk+l/y+8Zdrif86CNg+BRapp7PpshomoYc4bbPuKiw8fXWaTfC4ed0Xee391z8TXh1EXzwTHR2t3BHfuOTwhhWD2Pi2jv6ye8T/+KWuzfmNyaRPHj07a3UVldw2vTR0ZVv3QNv/gYGeLO8nfF/ChOciARC/+vwFmqDYZPhGjfs1/drb+TBSAOvfHkBDApgbcAZN3Vfd8wl8Oi1XSe/8JPfQz+ev9ikcAYOgy8vc/d3vAe3zo+eD7s3FCwskb7U2hFmycrtXDCrjqqK8uiGzgNQMx6+vqpwwYlIYPS/mt9Q28HmDm2dYR59eyvnzhzHsCAmvj2Jbc8ZO82tmj2UHr/zm98O3D8Xgti2XaQHS1fvYH9HmAtn13XdoCmNRaQX+l/y2xnt9PXUeztoaQvxmeP64fBfsR3glPyWNv8992v//Y5vvZk9TiQAHnl7K2NrqzhxysiuGzSlsYj0Qv9LfkPR4b4eeqOBuqHVzD90ZIonBVC3mt/27uulNFQkqfkV6Uf2HOjkr2saueCYOsrL4q5qhNpV9olI2vpX8tu6GzY8D+WVbN/TxvPvN/GpYyd2Lyj7g9gpjBvegDvP7r5eSoM/qsmSb8Fdn4BIp3s8ZHTy54gEzBMrttEZtt3H9n3wC7DuaZV9IpK2/pX87tniluOO4eG3thCx8OnjJhY2pr5y5r/Byde5+9vfhbY9cMI/wCGnFDYuyb/KgXD+j2H8bGh4Faq96bvP/s/CxiWSQ4+8vZWpowYzc0Lc5EVb3oDR0+H0GwsTmIgETv9Kfv3LvlNO5eG3Gjhu8nCmjBpc2Jj6yuFnw5nfhYHDobXZrZv/ZagaUti4pDDmfhGmLXRD3YU7YN41Ohek32jc28bL6z/i/Fl1mPiOnJEwTJgLU04tTHAiEjgpk19jzJ3GmEZjzIqYdSOMMU8ZY973lsP7Nsw0eW0dG/Za1u7Yx4Wz6lI8oR+oGAid+6P3pXT5vd07D6jne4kLVLmdhiUrd2AtnHf0uO4bIyEo73+jdopI30mn5vduIG6WBW4EllprDwOWeo8LL+R6u7+4cS/GwLkzExSU/U1sJw91+Chtsb3d1fO91N1NUMrtNCxesZ0powYzfWxN942REJQp+RWR9KVMfq21zwG74lZfBPzau/9r4OIcx5W+jz6ArW+5oW42vwrAsx/u5fhDRjCmtgQSACU84tMPIfEUfbndC7v3d7Dsw484d+a47k0eQMmviPRapm1+x1prtwF4yzHJdjTGXG2Med0Y83pTU1OGL9eD206BRafDG3fDs98D4M2dFZx/zPjcv1YxGuId+qqhSnhKXeyPn8FJP5JSutIqt/u8zO6lp1btIByxnDszSZkeCSv5FZFe6fMSw1q7CFgEMHfu3NyPut95wC1b3EgPdx77MDuWtXJOKTR5APjsPW462yFjNaOXOONnw6xLCx2FBFSfl9m99MTybUwcPrD7KA++cKeSXxHplUxrfncYY8YDeMvG3IWUodZmMGXcv67cNXmoKZEmANW1MP4YqBlb6EikWIw+Qj+EJJHiK7dTaGnr5IV1OznnqCRNHkDNHkSk1zJNfh8FrvDuXwE8kptwstC2h0hFNWsa93HWDCWCUsKU+EpixVdup/CXVY10hi3nJhrlAdwU3lbNHkSkd9IZ6ux3wDJgujGmwRhzFfBfwJnGmPeBM73HhdXWTAcDAFhwpJJfKWG24FeqpcACU26n8MSKbYytrWLOpCSjskXCbqmhzkSkF1KWGNbay5JsWpDjWHqvfW/0/rZ3ORCpYOrowf13YgsRkTQUdbmdpgMdIZ5d28Rn506iLNkU9ZGQW6rmV0R6IdgzvO147+Dd8OAxPNZxHAtV6yulavp5MGkenP7NQkcikrXn1jbR1hnh7J46L4c73LKsMj9BiUi/EOyfy96kFlz5BEtapvDte9/k90doiCcpUQOHwVVPFjoKkZx46r1Ghg6s5IRDRiTfyZvVUzMaikhvBLvmN9TulhVVPL3KFZRzJwdmxk4REUkgHLE8s6aRM6aPpqK8h68pP/nVBD8i0gvBTn47Xc2vrajm2bVNnHp4ioJSRESK3lubdrNrfwcLU43c06nkV0R6L9iZolfz+8HuEDv3tfOxw0YVOCAREcnWU6t2UFFmOPXw0cl3sjY6yZGSXxHphWC3+f1gKQAvb3a//pX8iogE39JVjZw4dQS11T10ZLvzbNj8irtfPiA/gYlIvxDsml+vwHtqs+HQ0YMZP3RggQMSEZFsbNi5n3WN+1KP3OMnvgDjZ/VtUCLSrwQ7+Q21YYdP4ZX1H/Gxw3q4PCYiIoHw9KodAL0btrJ6aB9FIyL9UeCT39ZIBW2dEU6ZpiYPIiJBt3RVI9PH1jBpxKD0n1RR1XcBiUi/E/Dkt52WUDnlZYYTp/YwFqSIiBS9PQc6eXXDLhYcmWK89nCo62OTZAY4EZEEgp38bnyJj9rLmTNpGDU9dYwQEZGi99e1jYQjNo0hzvbnJyAR6ZeCO9pD82bo2Ecksp+TDh1Z6GhERCRLT69qZNSQAcyeOKznHbe+Hb1fU9e3QYlIvxPc5Le9BYBfhc7nkilq8iAiEmSd4Qh/XdPIuTPHUVaWohmDN8ERl/8BJs/v++BEpF8JbrMHb2af/WYgx9ZrSmMRkSB7bf0u9raFWJDOKA/+tMa1dTBgcN8GJiL9TnCTX6/wqxs5jMFVwa3AFhER1+RhQEVZepMV+cmvRnkQkQwENvntaHfTWk6rU3tfEZEgs9by9KodnHzoSAYNSKMy42Dyq2mNRaT3Apv87ljjZveZPqkXA6GLiEjRWde4j027DvDJSQfgvs/Cg1+A9n3RHRpeh99fDrvWw4YX4eXb3PpKzeopIr0X2OS3sWknANNnzC5wJCIiko2nvFndTqtcCWsXw8qHoWl1dIfX7oDVj8GGF2DVo7BzLRz1SahOMSqEiEgCgW0s29Tcwn4GMny4RnoQEQmypasamTmhlqFVMfUxftMGgEgoui7UBoNHwSV35zVGEek/AlnzG4lYmlv2YtXZQUQk0Hbua+fNTbtZcMTYaJILXZPfg+va3U1tfUUkC1nV/BpjNgB7gTAQstbOzUVQqXzQtI/ycBtmkNp7iYj0RqHK7WSeWd2ItXDmjLGwPib57UyU/La6MX6V/IpIFnLR7OEMa+3OHPyd9LTvY+W6DdSYViqrlPyKiGQgv+V2D5auamRcbTVH1dXCupjkt63ZLa2FA16orc3QvldDnIlIVoLV5nf7cvjlaVxsw1AOdqA6u4mIBFV7KMzz7zdx0ZwJGGO6Nnt45Bo44nx4+Vb44C9u3bJb3HLyyfkPVkT6jWyTXwssMcZY4JfW2kXxOxhjrgauBqivr8/u1Zo3gw3zYNWnaB04js+ff0l2f09EpPT0WG7ntMxO4eUPd7G/I8zCI8e4FZEQmHI4/BxY82fYvxN2fei2LfjX6Gxu9ZrSWEQyl23ye7K1dqsxZgzwlDFmtbX2udgdvIJ1EcDcuXNtVq/mdYC4fe88zjrudJgwPas/JyJSgnost3NaZqewdNUOqivLOOlQb1a3SCeUVcCsS13y64/uMPoI+Nj1fRmKiJSQrEZ7sNZu9ZaNwMPACbkIKqlQOwAHbCWzJ2l8RxGR3sp7uZ08DpauauSUaaOprix3KyNhKK+MdmjT6A4i0gcyTn6NMYONMTX+feAsYEWuAkso1ApAux2g5FdEpJcKUm4nsXr7XrY0t3LmjDHRlZEQlJVDpZ/8trqbkl8RyaFsmj2MBR42xvh/5z5r7eKcRJXMppcBGD28lpFD1NtXRKSX8l9uJ7HUm9XtrMHvwytPuJXb3nXNHvxkd88WV/OraYxFJIcyTn6ttR8Cs3IYS2qbXwFg2qTxeX1ZEZH+oCDldhJPr2pk1qRhDP/zF2B/Y3RD3RyorXP3m1a7Nr8DhxckRhHpnwI11Fk4HOaJ8InMnDiy0KGIiEiGmva2805DM9cvPBxe3gfHfwnO+D9uY1WNa/eLARtxk11oXF8RyaFATW8c6mijxQ5yg6GLiEgg+bO6LThijFezOwwGjXC38kq3U1WN1+GtDSrU7EFEcidQya8NtdHOAI6qG1roUEREJENPrdpB3dBqjhxb7Wp3E3Voq6j2hjprV82viORUoJLf8nAbA6oHMXRQZaFDERGRDOxvD/Hc2ibOOmocxhu+sufkt1Ud3kQkp4KT/LbtodJ2MqxmSKEjERGRDD2zppH2UIRzZ46D5k1uZaKa3YoqWPsktLWo5ldEciowHd72b1vNYGDoMPX6FREJqidWbGfUkAHMPWQErHUj+DB4dPcdj78K3l8CpgymfyK/QYpIvxaY5HdDYwtHATX1swsdioiIZKCtM8wzqxu5eM4EysvMwYmLGH1E953n/ZO7iYjkWGCaPWxobAGgfnRNgSMREZFMPLe2iQMdYdfkAQ5OWa9mDSKST4FJfjfu3APA8CGDChyJiIhkYvGK7QwdWMm8qd5Y7aE2t9T0xSKSR4Fp9rB55153p0wjPYiIBE1b637eW7WcS6aNoXLfNoh0Qss2t7FSya+I5E8gkt+OUITG5n1QiZv3XUREAmXfnZ9iMS/DOuDHMRtMuSaxEJG8CkQm+eHOfRgbdg/KygsbjIiI9Fq4eTPLzeEcbde6Fef/xLX1HTpJNb8ikleBSH7XbN9LOX7yG4iQRUTEs3t/B5GOVsKj58BOL/md8zkoV3kuIvkXiA5va7bvpaos4h6Uq82viEiQPPbuVqrooH7MiOhKJb4iUiCBSH7X7tjL+Bov6VXNr4hIYFhrufeVTQw0IYbXaqhKESm8QGSSQ7a+wDmVr7kHavMrIhIYHz5zNwuaXqKqsgNTqY5tIlJ4RZ/87msPcV3rrUxp3wEDR8CgkYUOSURE0tFxgEOf+yrfqASLgdHT4cgLYO/2QkcmIiWs6JPf93fspc600zDlM0z83K+gLBAtNURESt6H23YyFXiq/mucecXNrp3v0ZcUOiwRKXFFn0mu3bGXajqorR2qxFdEJEDuenYVAPOPrI92cDPG3URECiSrbNIYc44xZo0xZp0x5sZcBRVr9fa9VNFJzeAhffHnRURKSj7KbYCX1u3khVWbARii8ltEikjGya8xphz4OXAuMAO4zBgzI1eB+dZub6HadGI0CLqISFbyVW5v3nWA6x94h6nDvNreiqpcv4SISMayafN7ArDOWvshgDHm98BFwHu5CMz33S1fcndUeIqIZKvPy+1Xf/q3DN/1DvcYqK/x6lcqVHkhIsUjm+R3ArA55nEDcGL8TsaYq4GrAerr63v1Ah2hCM2102moOZqJ0z+RRagiIkIa5XY2ZTZAuHYSHfYAh4wazIAB5VB5Ckw8PouQRURyK5vkN1GPBdtthbWLgEUAc+fO7ba9JwMqyjju+j9kFp2IiMRLWW5nU2YDzL/ye5lFJiKSJ9l0eGsAJsU8nghszS4cERHpQyq3RaTkZZP8vgYcZoyZYowZAFwKPJqbsEREpA+o3BaRkpdxswdrbcgYcy3wJFAO3GmtXZmzyEREJKdUbouIZDnDm7X2ceDxHMUiIiJ9TOW2iJQ6TZkmIiIiIiVDya+IiIiIlAwlvyIiIiJSMpT8ioiIiEjJMNb2egzzzF/MmCZgYwZPHQXszHE42VJM6Sm2mIotHlBM6Sp0TJOttaML+Pp5pzK7zymm9Cim9BRbTMUQT8JyO6/Jb6aMMa9ba+cWOo5Yiik9xRZTscUDiildxRiTJFaM75ViSo9iSo9iSq3Y4omlZg8iIiIiUjKU/IqIiIhIyQhK8ruo0AEkoJjSU2wxFVs8oJjSVYwxSWLF+F4ppvQopvQoptSKLZ6DAtHmV0REREQkF4JS8ysiIiIikjUlvyIiIiJSMoo6+TXGnGOMWWOMWWeMubFAMUwyxjxjjFlljFlpjLnOWz/CGPOUMeZ9bzm8ALGVG2PeMsY85j2eYox5xYvpfmPMgDzHM8wY85AxZrV3vOYX+jgZY77mvW8rjDG/M8ZU5/s4GWPuNMY0GmNWxKxLeFyM8z/eOf+uMebYPMb03957964x5mFjzLCYbTd5Ma0xxpydr5hitt1gjLHGmFHe47wcJ+m9QpfbKrN7FY/K7MQxqMzOMKaYbUVdZhdt8muMKQd+DpwLzAAuM8bMKEAoIeDr1tojgXnANV4cNwJLrbWHAUu9x/l2HbAq5vH3gB97Me0GrspzPD8FFltrjwBmebEV7DgZYyYAXwHmWmtnAuXApeT/ON0NnBO3LtlxORc4zLtdDdyax5ieAmZaa48B1gI3AXjn+6XAUd5zfuF9PvMRE8aYScCZwKaY1fk6TtILRVJuq8xOn8rsxO5GZXamMQWjzLbWFuUNmA88GfP4JuCmIojrEdybugYY760bD6zJcxwTcR/AjwOPAQY3k0pFouOXh3hqgfV4nShj1hfsOAETgM3ACKDCO05nF+I4AYcAK1IdF+CXwGWJ9uvrmOK2fRK417vf5bMHPAnMz1dMwEO4L+YNwKh8HyfdevX+FV25rTI7aTwqs3uORWV2hjEFocwu2ppfoh8CX4O3rmCMMYcAc4BXgLHW2m0A3nJMnsP5CfAvQMR7PBJottaGvMf5Pl5TgSbgLu+y3u3GmMEU8DhZa7cAP8D9+twG7AHeoLDHyZfsuBTLef9F4AnvfsFiMsZcCGyx1r4Tt6lYjpN0VVTvi8rsHqnM7h2V2WkISpldzMmvSbCuYOOyGWOGAH8AvmqtbSlUHF4s5wON1to3Ylcn2DWfx6sCOBa41Vo7B9hPYS4rHuS1yboImALUAYNxl17iFdN4f4V+HzHGfAt36fhef1WC3fo8JmPMIOBbwLcTbU6wrpjex1JVNO+LyuyUVGbnRqHfR5XZGSjm5LcBmBTzeCKwtRCBGGMqcYXovdbaP3qrdxhjxnvbxwONeQzpZOBCY8wG4Pe4y2g/AYYZYyq8ffJ9vBqABmvtK97jh3AFayGP00JgvbW2yVrbCfwROInCHidfsuNS0PPeGHMFcD5wufWuTRUwpkNxX4LveOf6ROBNY8y4AsYkPSuK90VldlpUZveOyuzUAlNmF3Py+xpwmNfLcwCu8faj+Q7CGGOAO4BV1tofxWx6FLjCu38Frl1ZXlhrb7LWTrTWHoI7Ln+x1l4OPAN8pkAxbQc2G2Ome6sWAO9RwOOEu3Q2zxgzyHsf/ZgKdpxiJDsujwKf93rGzgP2+Jfa+pox5hzgm8CF1toDcbFeaoypMsZMwXVYeLWv47HWLrfWjrHWHuKd6w3Asd65VrDjJD0qeLmtMjvtmFRm947K7BQCVWYXqrFxOjfgPFwPxg+AbxUohlNwVfPvAm97t/Nw7bWWAu97yxEFiu904DHv/lTcCb4OeBCoynMss4HXvWP1J2B4oY8T8G/AamAF8FugKt/HCfgdrv1aJ64wuCrZccFdGvq5d84vx/V6zldM63Btsvzz/LaY/b/lxbQGODdfMcVt30C080RejpNuGb2PBS23VWb3KhaV2YljUJmdYUxx24u2zNb0xiIiIiJSMoq52YOIiIiISE4p+RURERGRkqHkV0RERERKhpJfERERESkZSn5FREREpGRUpN5FpDCMMf7QMgDjgDBuOk6AA9bakwoSmIiIdKMyW4JCQ51JIBhjvgPss9b+oNCxiIhIz1RmSzFTswcJJGPMPm95ujHmWWPMA8aYtcaY/zLGXG6MedUYs9wYc6i332hjzB+MMa95t5ML+x+IiJQOldlSTJT8Sn8wC7gOOBr4HHC4tfYE4Hbgn719fgr82Fp7PPBpb5uIiOSfymwpKLX5lf7gNevNEW6M+QBY4q1fDpzh3V8IzHDTxQNQa4ypsdbuzWukIiKiMlsKSsmv9AftMfcjMY8jRM/xMmC+tbY1n4GJiEg3KrOloNTsQUrFEuBa/4ExZnYBYxERkZ6pzJY+o+RXSsVXgLnGmHeNMe8B/1jogEREJCmV2dJnNNSZiIiIiJQM1fyKiIiISMlQ8isiIiIiJUPJr4iIiIiUDCW/IiIiIlIylPyKiIiISMlQ8isiIiIiJUPJr4iIiIiUjP8PqDVqnCKodToAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "#Plot the results\n",
    "plt.figure(figsize = (12, 4))\n",
    "plt.subplot(121)\n",
    "plt.title(\"Transcript T\")\n",
    "plt.plot(timepoints, results_det[\"T\"], label = \"deterministic\")\n",
    "plt.plot(timepoints, results_stoch[\"T\"], label = \"stochastic\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Time\")\n",
    "\n",
    "plt.subplot(122)\n",
    "plt.title(\"Protein X\")\n",
    "plt.plot(timepoints, results_det[\"X\"], label = \"deterministic\")\n",
    "plt.plot(timepoints, results_stoch[\"X\"], label = \"stochastic\")\n",
    "plt.legend()\n",
    "plt.xlabel(\"Time\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### More examples / assignments:\n",
    "\n",
    "1. Can you obtain a reduced model with Hill functions for the signaling system discussed earlier? Using the various techniques discussed in this lecture, write that model using bioscrape. Simulate and compare with full model dynamics (with massaction propensities) from Lecture 1 [here](http://www.cds.caltech.edu/%7Emurray/courses/be240/sp2020/W2_bioscrape.ipynb).\n",
    "\n",
    "2. For the repressilator model in [[1](https://www.nature.com/articles/35002125)], write a bioscrape model using non-massaction propensities\n",
    "\n",
    "3. A transcription-translation model using non-massaction propensities is available in `Basic Examples --- START HERE.ipynb` in the `bioscrape/examples` directory. \n",
    "\n",
    "4. Another such model is given next. Feel free to explore this code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Bioscrape model #rxns= 2 \n",
      "rxns= [(['D', 'A'], ['E'], 'hillpositive', {'k': 0.295, 'K': 0.69, 'n': 2.166, 's1': 'F'}), (['E', 'F'], ['D', 'G'], 'hillpositive', {'k': 1.19, 'K': 1.24, 'n': 0.143, 's1': 'E'})]\n",
      "Warning! SpeciesB not currently used in any rules or reactions.\n",
      "Warning! SpeciesC not currently used in any rules or reactions.\n",
      "Loaded Model \n",
      " [(<bioscrape.types.PositiveHillPropensity object at 0x0000020B4A0CE990>, <bioscrape.types.NoDelay object at 0x0000020B4A124170>, {'D': -1, 'A': -1, 'E': 1}, {}), (<bioscrape.types.PositiveHillPropensity object at 0x0000020B4A10EA20>, <bioscrape.types.NoDelay object at 0x0000020B4A1240B0>, {'E': -1, 'F': -1, 'D': 1, 'G': 1}, {})]\n"
     ]
    }
   ],
   "source": [
    "## Example not shown in lecture\n",
    "\n",
    "import numpy as np\n",
    "import pylab as plt\n",
    "from bioscrape.simulator import *\n",
    "from bioscrape.types import *\n",
    "import warnings\n",
    "\n",
    "\n",
    "#Parameter ranges to randomly choose parameters (on a log scale)\n",
    "param_min = -2\n",
    "param_max = 2\n",
    "\n",
    "#Names of different supported propensities\n",
    "propensity_types = ['massaction', 'hillpositive', 'proportionalhillpositive', 'hillnegative', \n",
    "                    'proportionalhillnegative', 'massaction']#, 'general']\n",
    "\n",
    "#parameter names required for each propensity (general will be treated by itself)\n",
    "propensity_param_requirements = {\n",
    "\t'massaction':['k'], 'hillpositive':['k', 'K', 'n'], 'hillnegative':['k', 'K', 'n'],\n",
    "\t'proportionalhillpositive':[\"k\", \"K\", \"n\"], 'proportionalhillnegative':[\"k\", \"K\", \"n\"]\n",
    "}\n",
    "#species (passed in as parameters) required for each propensity (general will be treated by itself)\n",
    "propensity_specie_requirements = {\n",
    "\t'hillpositive':['s1'], 'hillnegative':['s1'], 'proportionalhillpositive':['s1', 'd'], \n",
    "    'proportionalhillnegative':['s1', 'd'], \"massaction\":[]\n",
    "}\n",
    "\n",
    "\n",
    "species = ['A', 'B', 'C', 'D', 'E', 'F', 'G']\n",
    "n_species = len(species)\n",
    "n_reactions = np.random.randint(1, 3)\n",
    "\n",
    "reactions = []\n",
    "for r in range(n_reactions):\n",
    "\n",
    "    try_again = True\n",
    "    while try_again:#Loop ensures no positive feedback which leads to long simulations\n",
    "        inputs = []\n",
    "        outputs = []\n",
    "        while(len(inputs) == 0 and len(outputs) == 0):\n",
    "\n",
    "            n_inputs = np.random.randint(0, 5)\n",
    "            for i in range(n_inputs):\n",
    "                inputs.append(species[np.random.randint(0, len(species))])\n",
    "\n",
    "            n_outputs = np.random.randint(0, 5)\n",
    "            for i in range(n_outputs):\n",
    "                outputs.append(species[np.random.randint(0, len(species))])\n",
    "\n",
    "        inputs_in_outputs = len([i for i in inputs if i in outputs])\n",
    "        if inputs_in_outputs >= len(inputs):\n",
    "            try_again = True\n",
    "        else:\n",
    "            try_again = False\n",
    "\n",
    "    prop_type = propensity_types[np.random.randint(0, len(propensity_types))]\n",
    "    param_dict = {}\n",
    "    if prop_type != 'general':\n",
    "        required_params = propensity_param_requirements[prop_type]\n",
    "        required_species = propensity_specie_requirements[prop_type]\n",
    "        param_dict = {}\n",
    "        for p in required_params:\n",
    "            param_dict[p] = round(np.exp(np.random.uniform(low = param_min, high = param_max)), 3)\n",
    "        for i in range(len(required_species)):\n",
    "            k = required_species[i]\n",
    "            param_dict[k] = species[np.random.randint(0, len(species))]\n",
    "\n",
    "    elif prop_type == 'general': #Here we will use a random(ish) rational function\n",
    "        rate_str = \"(1+\"\n",
    "        numerator_terms = np.random.randint(0, 5)\n",
    "        denominator_terms = np.random.randint(0, 5)\n",
    "        for i in range(numerator_terms):\n",
    "            coef = str(round(np.exp(np.random.uniform(low = param_min, high = param_max)), 3))\n",
    "            exp = str(round(np.random.uniform(low = 0, high = param_max), 3))\n",
    "            specie = species[np.random.randint(0, len(species))]\n",
    "            rate_str += coef+\"*\"+specie+\"^\"+exp+\"+\"\n",
    "        rate_str =  rate_str[:-1] + \")\"\n",
    "        rate_str += \"/(1+\"\n",
    "        for i in range(denominator_terms):\n",
    "            coef =str(round(np.exp(np.random.uniform(low = param_min, high = param_max)), 3))\n",
    "            exp = str(round(np.random.uniform(low = 0, high = param_max), 3))\n",
    "            specie = species[np.random.randint(0, len(species))]\n",
    "            rate_str += coef+\"*\"+specie+\"^\"+exp+\"+\"\n",
    "        rate_str =  rate_str[:-1] + \")\"\n",
    "        param_dict['rate'] = rate_str\n",
    "\n",
    "    rxn = (inputs, outputs, prop_type, param_dict)\n",
    "    reactions.append(rxn)\n",
    "\n",
    "\n",
    "print(\"Bioscrape model #rxns=\", len(reactions), \"\\nrxns=\", reactions)\n",
    "M = Model(reactions = reactions, initial_condition_dict = {s:np.random.randint(10, 100) for s in species})\n",
    "M.write_bioscrape_xml('lecture2_model.xml')\n",
    "M2 = Model('lecture2_model.xml')\n",
    "print(\"Loaded Model \\n\", M2.get_reactions())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Rules in bioscrape\n",
    "### Adding Rules\n",
    "In deterministic and stochastic simulations, bioscrape also supports rules which can be used to set species or parameter values during the simulation. Rules are updated every simulation timepoint - and therefore the model may be sensitive to the timepoints spacing.\n",
    "\n",
    "For example:\n",
    "\n",
    "1. $I = I_0 H(T)$ where $H$ is the step function. Represents the addition of the inducer I at concentrion $I_0$ some time T. Prior to t=T, I is not present.\n",
    "2. $S = M \\frac{X}{1+aX}$ represents a saturating signal detected from the species X via some sort of sensor.\n",
    "\n",
    "To implement these rules in bioscrape:\n",
    "```\n",
    "# Usage: rule_tuple = (rule_type, rule_dict, rule_frequency)\n",
    "\n",
    "rule1 = (\"assignment\", {\"equation\": \"I =_I0*Heaviside(t-_T)\"}, \"repeated\")\n",
    "\n",
    "rule2 = (\"assignment\", {\"equation\": \"S = 50*X/(1+.2*X)\"}, \"repeated\")\n",
    "\n",
    "M = Model(rules = [rule1, rule2])\n",
    "\n",
    "# Or use create_rule function\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "Rules can also be used for quasi-steady-state or quasi-equilibrium approximations, computing parameters on during the simulation, and much more!\n",
    "\n",
    "There are two main types of rules:\n",
    "1. \"additive\": used for calculating the total of many species. Rule 'equation' must be in the form $s_0 = s_1 + s_2 ...$ where $s_i$ each represents a species string.\n",
    "2. \"assignment\": a general rule type with 'equation' of the form $v = f(s, p)$ where $v$ can be either a species or parameter which is assigned the value $f(s, p)$ where $s$ are all the species and $p$ are all the parameters in the model and $f$ is written as an string.\n",
    "\n",
    "\n",
    "For more information on Rules in bioscrape, refer to the [bioscrape Wiki](https://github.com/ananswam/bioscrape/wiki) and the `Basic Examples --- START HERE.ipynb` notebook in the `bioscrape/examples` directory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "source": [
    "# Bonus content: Local sensitivity analysis\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "slideshow": {
     "slide_type": "skip"
    }
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "kernelspec": {
   "display_name": "Python 3.7.4 64-bit",
   "language": "python",
   "name": "python37464bita6c5c53d7b16474cae5bc645c5b09138"
  },
  "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.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}