docs/source/examples/lecture2_michaelis_menten.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "-XqS9Lv7LRYS"
   },
   "source": [
    "# Michaelis Menten"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "cnn37173LU83"
   },
   "source": [
    "The Michaelis-Menten (MM) kinetic describes the rate of a chemical reaction catalized by enzymes, varying its concentration and the concentration of the substrate.\n",
    "\n",
    "In biochemistry, Michaelis–Menten kinetics is one of the best-known models of enzyme kinetics.\n",
    "\n",
    "It is named after German biochemist Leonor Michaelis and Canadian physician Maud Menten.\n",
    "\n",
    ">**Def:** Catalysis is the process of increasing the rate of a chemical reaction by adding a substance known as a catalyst.\n",
    "Catalysts are not consumed in the catalyzed reaction but can act repeatedly.\n",
    "\n",
    ">**Def:** Enzymes are proteins that act as biological catalysts (biocatalysts). \n",
    "Catalysts accelerate chemical reactions. \n",
    "The molecules upon which enzymes may act are called substrates, and the enzyme converts the substrates into different molecules known as products. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RH_7Kl56LWCb"
   },
   "source": [
    "The MM kinetic is model as\n",
    "\n",
    "$$\n",
    "E + S \\underset{k_b}{\\stackrel{k_f}{\\rightleftharpoons}} ES \\xrightarrow[]{k_c} E + P\n",
    "$$\n",
    "\n",
    "where $k_c$ is the **catalytic constant** which regulates the rate of the enzymatic reaction.\n",
    "\n",
    "The link between enzyme and substrate allows to reach the activation energy of the (*irreversible*) reaction.\n",
    "\n",
    "The enzyme is highly specific for its substrate and their link is called **lock and key** mechanism.\n",
    "This mechanism allows the union between enzymes and substrates non properly associated, leading to a consequent loss in the reaction efficiency.\n",
    "In this way we can associate to each enzyme an affinity distribution with its substrates."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "rEN24GtRLZOP"
   },
   "source": [
    "Using the **action mass law** we can re-write the MM kinetic model with the following set of differential equations\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "    \\frac{d[E]}{dt}  = -k_f[E][S] + k_b[ES] + k_c[ES] \\\\\n",
    "    \\frac{d[S]}{dt}  = -k_f[E][S] + k_b[ES]           \\\\\n",
    "    \\frac{d[ES]}{dt} =  k_f[E][S] - k_b[ES] - k_c[ES] \\\\\n",
    "    \\frac{d[P]}{dt}  =  k_c[ES]                       \n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    ">⚠️**Note:** The 2nd equation is exactly the opposite of the 1st one.\n",
    "> This imposes the conservation of enzyme concentration!\n",
    "\n",
    ">⚠️**Note:** The 4th equation is a typical 1st order irreversible reaction."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ClFqAtMgLbF7"
   },
   "source": [
    "The solution of this system of equations can be obtained using the *Stationary State approximation*, i.e\n",
    "\n",
    "* The 1st equation is redundant, since the kinetic imposes the conservation of the enzyme concentration\n",
    "\n",
    "$$\n",
    "[E]_{tot} = [E] + [ES] \\rightarrow [E] = [E]_{tot} - [ES]\n",
    "$$\n",
    "\n",
    "* The approximation of (pseudo) Stationary State assumes as constant the concetration of $[ES]$, i.e more is created and more is consumed\n",
    "\n",
    "$$\n",
    "\\frac{d[ES]}{dt} = 0 \\rightarrow [ES] = \\frac{k_f}{k_b + k_c}[E][S]\n",
    "$$\n",
    "\n",
    "* Therefore, we can define a reaction constant given by the combination of $k_f$, $k_b$, and $k_c$, i.e\n",
    "\n",
    "$$\n",
    "k_m = \\frac{k_b + b_c}{kf}\n",
    "$$\n",
    "\n",
    "called **Michaelis-Menten constant**.\n",
    "\n",
    "By substitution in the previous formula, we obtain\n",
    "\n",
    "$$\n",
    "[ES] = \\frac{1}{k_m}[S]([E]_{tot} - [ES]) = \\frac{[E]_{tot} [S]}{k_m + [S]}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QQ5h5cFOLd6Q"
   },
   "source": [
    ">⚠️**Note:** The approximation of stationary state is considered *pseudo* since in the final result remains the dependency by $[S]$.\n",
    "> This approximation can be rethought as a very slow reaction process."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "jjHDDqNkLgGI"
   },
   "source": [
    "Now we can define the velocity of the product reaction as\n",
    "\n",
    "$$\n",
    "v_p = \\frac{d[P]}{dt} = - \\frac{d[S]}{dt} = \\frac{k_c [E]_{tot} [S]}{k_m + [S]}\n",
    "$$\n",
    "\n",
    "since in an *ideal* kinetic the substrate is consumed as much as the product is created.\n",
    "\n",
    "It is interesting to notice that *if* $[S] = k_m$ we are at the half of the maximum velocity of the reaction: in this way we can think about $k_m$ as the *affinity degree* between substrate and enzyme.\n",
    "\n",
    "If $[S]\\rightarrow\\infty$ we achieve the asymptotic value of the reaction velocity: it corresponds to the condition in which there are more substrate than necessary and therefore the reaction is as fast as possible.\n",
    "\n",
    "$$\n",
    "v_p = \\frac{v_{MAX} [S]}{k_m + [S]}\n",
    "$$\n",
    "\n",
    "where $v_{MAX}$ is the maximum velocity of the product generation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "FDmmGi-iLimw"
   },
   "source": [
    "Now lets try to implement the MM model as Python code.\n",
    "As before we can easily obtain a simulation of our model using the *odeint* function of *scipy*."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "executionInfo": {
     "elapsed": 343,
     "status": "ok",
     "timestamp": 1632140953776,
     "user": {
      "displayName": "nico nico",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "04264570794863623616"
     },
     "user_tz": -120
    },
    "id": "pUW2WrlgLSNU"
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def conversion (y : tuple, t : np.ndarray, \n",
    "                kf : float, kb : float, kc : float) -> tuple:\n",
    "    '''\n",
    "    Michaelis-Menten conversion\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "        y : tuple\n",
    "            List of [S], [E], [ES], and [P] arrays\n",
    "            \n",
    "        t : array-like\n",
    "            List of time points\n",
    "        \n",
    "        kf : float\n",
    "            Forward constant\n",
    "        \n",
    "        kb : float\n",
    "            Backward constant\n",
    "        \n",
    "        kc : float\n",
    "            Product constant\n",
    "    \n",
    "    Returns\n",
    "    -------\n",
    "        delta : tuple\n",
    "            Next point of [S], [E], [ES], and [P] arrays\n",
    "    '''\n",
    "    \n",
    "    S, E, ES, P = y\n",
    "    d_S  = -kf * S * E + kb * ES\n",
    "    d_E  = -kf * S * E + kb * ES + kc * ES\n",
    "    d_ES =  kf * S * E - kb * ES - kc * ES\n",
    "    d_P  = kc * ES\n",
    "    return (d_S, d_E, d_ES, d_P)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 420
    },
    "executionInfo": {
     "elapsed": 1071,
     "status": "ok",
     "timestamp": 1632140990773,
     "user": {
      "displayName": "nico nico",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "04264570794863623616"
     },
     "user_tz": -120
    },
    "id": "Yy1Lp9x9LkgB",
    "outputId": "23915131-a991-4822-e972-bee7a3c7abd9"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7eff4c7b3fd0>"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAfwAAAGBCAYAAACD7gMMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeZyNdf/H8df3zL4xY4YxjC3ZyZRJdacium/RYqksGWVJJUJKac+SKBVZCiGyRaXiViRU6k5ExlIq0swwxljG7Ms5398f1yE/WWauOWeuM2c+z8fjPM7MOd/r+r677wcf13V9F6W1RgghhBDezWZ1ACGEEEK4nxR8IYQQogKQgi+EEEJUAFLwhRBCiApACr4QQghRAUjBF0IIISoAX6sDuFNUVJSuW7eu1TGEEEKIMrFt27Z0rXXV833n1QW/bt26bN261eoYQgghRJlQSh280HdyS18IIYSoAKTgCyGEEBWAFHwhhBCiApCCL4QQQlQAUvCFEEKICkAKvhBCCFEBSMEXQgghKgAp+EIIIUQFIAVfCCGEqAA8ruArpeYqpdKUUrvO+qyKUmqdUuo353uElRmFEEKI8sbjCj4wH+h4zmdPAeu11g2A9c7fhRBCCFFMHreWvtb6a6VU3XM+vhNo6/z5PWAj8GSZhQJI+wWyUuGytpdqKYQQwotohwPsdrTWYLeDw2F8dta78bMG7fxda6P96Z8dDtAY32t95nufkBD8atYsk/8Ojyv4FxCttT7s/DkViL5QQ6XUIGAQQO3atV2X4JvJ8Nf3MGLXpdsKIYQoEa01Oj8fR04OjpxcdG4OjtxcHLl56Dzj3ZGXi84vMNrl5xk/Fxi/68ICHAXO3wsLz7woLEQXFKKLis56FUKR3fjZbofT73Y72vnirHd3Cu3QnlrTprm1j9PKS8E/Q2utlVL6It/PAmYBxMfHX7BdiUU1gMQPoCAH/INddlohhPAWurCQohMnsJ84if3ECewZGdgzTmLPyMCRkYH9VCaOrEzjPTMTe3YWjuwcHFlZOLKzjavhkvLzw+bnh/L3//t1+ndfX+NnPz9swUHg64vy9TM+9/UFXx+Ujy/K1wd8fFA2H5SfL9h8UD42493X5+/flQ18bCibD9hsKJsCmw/YFMrHB1DG90oZbZUyvrM5f8b5u1LG78qGb3Q1V//fcEHlpeAfUUrFaK0PK6VigLQyTxB5ufF+/A+o3qLMuxdCCKvYs7IpOnyIwrQ0itKOUpSWZrzS0yk6lo49/RhFx47hyMy88En8/PCpVAmf0FBslSrhExaKb7Vq2EJDsYWEGK/gYOM9KAhbSDAqMBBbUDC2oEBUQCC2wADjswDjXfn7G8VUFEt5KfifAvcBrzjfPynzBFENjPf036TgCyG8iiMvj4KDf1GYnERhcjIFSckUJidTeOgQhampOE6d+scxtkqV8I2KwjcykoAmjQmJjMI3sgo+ERH4hEfgEx6OT0Q4PpUr41OpEiooyLiyFZbxuIKvlFqCMUAvSimVDLyAUeg/UEoNAA4C95R5sCr1jfdjv5d510II4QpFx4+T/9vv5P/xOwW//0HBnwfI//NPig6nGgPJnGzBwfjFxuJXsybB8fH41YjBt3oMftWj8a1WDd+qVbEFBVn4XyLM8LiCr7XudYGv2pdpkHP5B0Pl2sYVvhBCeDBtt1Owfz95v/xC3t5fyP9lL3m//Ir9+PEzbWyhofjXq0dwq3j869bBv05d/GvXwq9WLXzCw+Vq3At5XMH3aFGXQ/o+q1MIIcT/U3joELk7dpC7M5HcXYnk7dmLzskBQPn7E9CgAaHt2hLQoAEBlzcgoMHl+FarJkW9gpGCXxKRDSBpkXHrS/6gCCEsoB0O8n/7nZwffyT3p5/I+eknilJTAaO4BzZpQni3bgS1aE5AkyYE1KuH8vOzOLXwBFLwSyKqARRkQWYqVIqxOo0QooIoSEoie/N3ZP/wP3J+2HLm1rxvdDTBra4i6MqrCLrqSgIbNpTiLi5ICn5JnJ6ad+w3KfhCCLdx5OWRs2ULWd98S/bXX1Nw8CBgFPjQG9oQfM21BLdujV/NGnJbXhSbFPySOHtqXr0brc0ihPAq9pMnydq0icwv15P17bfo3FxUQADB17Qm4t57CWnTBv96daXAC9Ok4JdEWA3wC5apeUIIlyg6cYLMdevI/Pxzsn/YAnY7vtWqUbnLnYS1a0dw69bYAgOtjim8hBT8krDZILK+TM0TQpjmyM4m88svyVi1muzvv4eiIvzq1CZywADCbulAYLNmsnqccAsp+CUV2QBStlmdQghRjmiHg5wtW8hY+Qmn1q5F5+TgV6MGkfffR9ittxLYtKncqhduJwW/pKIawO6PoTAP/ORWmxDiwoqOHuXkhx9xcvlyClNSsIWEUKnTrYR36UJQq1ZS5EWZkoJfUlENAQ3H90N0U6vTCCE8jNaanO+/58SSpWRu2ABFRQRfcw1Vhw8nrEN7WZJWWEYKfkmdPTVPCr4QwsmRk0PGp59x/P2FFPz+Bz4REVTp25fwu+8ioF49q+MJIQW/xE4XfBm4J4QACtPSOLFwISc+WI4jI4OApk2ImTCBSp07YfP3tzqeEGdIwS+pgFBjep5MzROiQss/cIDjc+eSsfITtN1O2C23UKVvAkFXXSXP5oVHkoJvRtTlcoUvRAWV9+uvpM+YSebatSg/Pyp370Zk//74165tdTQhLkoKvhmRDSBxhWyiI0QFkvfLL6RPn0HmunXYQkOJHDSIKgl98I2KsjqaEMUiBd+MqAaQnwHZRyG0mtVphBBulP/bbxydOpXMdV9iCw0lavBgqtzXF5/Kla2OJkSJSME3I/KsNfWl4AvhlQpTUjg6bToZn3yCLSiIqEceoUrfBCn0otySgm9G1FlT8+peb20WIYRL2U+eJH3mTE4sXgJKUeW++4gc9AC+ERFWRxOiVKTgm1G5FvgEyMA9IbyILijgxJIlHJ0xE0dmJpW7dqHqkCH4xchW2MI7SME3w+ZjzMeXqXlClHtaa7LWr+fIq69SePAvQq6/nmpPjiKwYUOrownhUlLwzYq6HA7vtDqFEKIU8vcf4Mi4cWR/9x3+9etTa9Y7hNxwg8yjF15JCr5ZdW+APZ/A/o1wWVuLwwghSsKRnU36229zbP572AICiH76aSJ69UT5+VkdTQi3kU2XzboyASrVhPVjjfn4Qohy4dS6dfzRqTPHZs+hcufO1P98DVX6JkixF15PCr5ZfoFw0yhI2Qr7Prc6jRDiEgoPHybpkSGkDH0Un/Bw6ixeTI1XJsjCOaLCkIJfGnH3QkQ9+GocOBxWpxFCnIe22zm+YAH7O99G9ubNVHviceqtWE7wVVdaHU2IMiUFvzR8/KDd03BkF+xZaXUaIcQ58vfv5+C9fTjy8gSCWrXislWfETlggNy+FxWSFPzSat4dqjaGDS+DvcjqNEIIQBcVcWzOHA506Ur+gQPUmDSRWrPewT821upoQlhGCn5p2Xyg3TPGqns7l1mdRogKL3//fv7sfS9pr00m9KYbqb/qMyrfcYdMtRMVnkzLc4Umt0NMHGycAE3vgIAwqxMJUeFoh4MT779P2uTXsQUFUfON1wnr2FEKvRBOcoXvCkpBx1fgVAr8d5TVaYSocAoPH+avAQM48vIEQq69lss++5RKt94qxV6Is0jBd5U618GNT8DPi2HncqvTCFFhZKxezf477iT3551UH/MSsW/PxLdqVatjCeFx5Ja+K904CvZvglUjIDYeqtSzOpEQXsuelc2R8ePJ+PhjguLiqPHqJPxr1bI6lhAeS67wXcnHF7rPBmWDDweCvdDqREJ4pdzEXRzo3o2MTz4havDD1Hl/oRR7IS5BCr6rhdeGO6YYK/BteNnqNEJ4Fa01x+bN589evdD5BdR5bz5VH30U5Ss3K4W4FPlT4g7NusIfX8G3b0Dd6+HyDlYnEqLcs588yaHRT5O1YQOhHdpTY9w4fMLDrY4lRLkhV/ju0nEiVGsKHz4AGclWpxGiXMvdsYP93bqR9e23RD89mti33pJiL0QJScF3F/9guGeB8Rx/+f1QVGB1IiHKHa01xxcs4M8+CShlo+6i96nSt69MtxPCBCn47hR1Odw5DZJ/hHXPW51GiHLFkZ3NoZEjOfLyBEJvvJF6H31I0BVXWB1LiHJLnuG7W7MukDQY/jcDal9jPN8XQlxU/v4DJD86lIL9B6g6YgSRDwxE2eT6RIjSkIJfFjq8BMlb4ZMhEN0cohpYnUgIj5X55ZccevIplL8/td+dQ8h111kdSQivIP9kLgu+/nD3fPANgGV9ID/L6kRCeBztcHB06lskDxmKf/361PvoQyn2QriQFPyyUrkm3DUX0vfBZ8NAa6sTCeEx7FlZJA8ZSvqMGVTu1o06CxfgFxNjdSwhvIoU/LJ0WVtjK91dK+DHOVanEcIj5B84wJ/39CBr0yain3mGmPHjsAUEWB1LCK8jz/DLWpvHjFH7n482ttStdbXViYSwTPZ335E8fATKZqP2u+8Scu01VkcSwmvJFX5Zs9mg69tQqQZ80Bey061OJIQlji9ezF8PDMIvOpq6K5ZLsRfCzaTgWyEoAnoshJxjsKI/OOxWJxKizOjCQlLHjOHImLGE3nADdZYsxj821upYQng9KfhWiWkJnSfDgU3w1Tir0whRJuyZmSQ9+BAnFi+hyoD+xE6fhk9oqNWxhKgQ5Bm+la5KgOQt8O3rEBsPjTtbnUgItylITiHpoQcp+PMgMePHE969m9WRhKhQ5Arfare+agze+/ghOPaH1WmEcIvcnTv5s0cPitKOUnvObCn2QlhACr7V/AKNTXZsPrAsAQpyrE4khEudWruWgwl9sQUFUXfpEkKuvdbqSEJUSOWq4CulRiildiuldimlliilAq3O5BIRdaDbHEjbA6tGyKI8wmscX7CAlGHDCWzcmLofLCPgssusjiREhVVuCr5SqibwKBCvtW4O+AA9rU3lQg06QNunYOdS2Pqu1WmEKBXtcHDklYkceXkCYR3aU/u9+fhWqWJ1LCEqtHJT8J18gSCllC8QDByyOI9r3TgKGvwb1jwFST9anUYIUxz5+aSMHMnx+fOJ6NOHmm++iS3QO27GCVGelZuCr7VOAV4D/gIOAxla67XWpnIxmw26vgOVYmD5fbIojyh37KdOkTRgIJlrPqfaqFFEP/M0ysfH6lhCCMpRwVdKRQB3AvWAGkCIUqrPedoNUkptVUptPXr0aFnHLL3gKtDjfeeiPP3AXmR1IiGKpfBIGgf7JJDz88/UmPwakf37oZSyOpYQwqncFHygA3BAa31Ua10IfAT869xGWutZWut4rXV81apVyzykS8S0hM6vw4GvYYMsyiM8X/7+Axzs1YvC5GRqv/M2lTvLmhJCeJryVPD/Aq5VSgUr47KhPbDX4kzuc+W90KoffPsG7P3M6jRCXFDuzz9zsHdvHPn51F64gJB//ePf4UIID1BuCr7W+gdgBfATkIiRfZalodzt1olQ4yr4+GFI/93qNEL8Q9bmzRzs1x9bWBh1Fy8iqFkzqyMJIS6g3BR8AK31C1rrxlrr5lrrBK11vtWZ3Mo3wFiUx9cflvWB/CyrEwlxxqnPvyDpoYfxr1WLuosX4V+njtWRhBAXUa4KfoUUXgvumgvpv8KnQ2VRHuERTnzwASmPPUZQixbUWfAevuV1vIwQFYgU/PLgsrZw83Ow+yP430yr04gK7ticOaQ+/wIhba6n9rtz8Klc2epIQohikIJfXrQZAY1vg3XPwcHvrE4jKiCtNWmvv0Haa5Op1KkTtaZNwxYUZHUsIUQxScEvL5SCLjMgoi4svx8yU61OJCoQ7XBwZPzLHJs1i/C776bGq5NQ/v5WxxJClIBvaQ5WSnUE7gKaAxGA30Waa611/dL0V+EFVjYW5Zl9M3xwH9z3mTGgTwg30nY7h599joyPP6ZKv35UG/WELKgjRDlkquArpfyAZRgr3wEU50+/jDZzhWpN4M5psKI/rH0WOk2yOpHwYrqwkJQnRpH5+edEDRlC1CODpdgLUU6ZvcJ/EuiCUcRXAyuBFCDPRbnExTTvDsnb4H/TITYerrjH6kTCCzkKCkgZPoKsr76i2qhRRPbvZ3UkIUQpmC3492IU+9Faa7nEtMItL8HhHfDpo8ZVf/UWVicSXsSRl0fykKFkf/st1V94nohevayOJIQoJbOD9uoCDuAt10URJeLjB3fNg6BwY1Ge3BNWJxJewpGdTdKDD5G9eTMx48dJsRfCS5gt+CeBTK11rivDiBIKizZW4stIgY8GgcNhdSJRztmzsvjrgUHkbN1KjUmTCO/e3epIQggXMVvwNwGVlVK1XBlGmFCrNXScAL+tha/l6Yowz56ZSdKAgeTu3EnNyZOpfPttVkcSQriQ2YI/DmOA3kQXZhFmXT0QWvaCja/Avi+sTiPKIfupU/w1YCC5e/YQ++YbVOr4H6sjCSFczFTB11rvwhil31EptUYp1VYpFeLaaKLYlILb3jAG7n30ABz7w+pEohyxnzzJX/36k793L7FTphDWoYPVkYQQbmCq4Cul7MDnQGXg38B64JRSyn6RV5ELc4tz+QVBj4WgbLAsAQqyrU4kygH7yZMc7N+f/N9+I3baW4Td3M7qSEIINzF7S1+ZeMkyvu4WURe6vwtpe4zperKznriI08W+4Pc/iJ0+jdCbbrI6khDCjczOw6/n0hTCdS5vD+2fg/VjoGYruG6w1YmEB7KfPMlf/Qf8XexvuMHqSEIINzNV8LXWB10dRLhQm8cg5Sdj6d2YK6BuG6sTCQ9yutjn//abFHshKhC5ze6NlIIuM6HKZcbOehkpVicSHsKekcFfAwaeeWYfeuONVkcSQpQRlxR8pVQ1pVRHpVSC89VRKVXNFecWJgVWgp6LoDAXPugLRflWJxIWs2dm8tfAB8jbt4+ab02VZ/ZCVDClKvhKqTZKqY3AYYxNdOY7X6uBw0qpDUqp60uZUZhVtZFxpZ+yFdaMsjqNsJA9K5ukBwaRt3cvsVPeJKxtW6sjCSHKmOmCr5R6CNgA3IAxCt8BpDlfdudnNwEblVIPlj6qMKXpHcYz/W3zYdt7VqcRFjDWxn+Q3MREak6eTNjNN1sdSQhhAbPz8K8EpgE+wGbgP0Co1jpGax0DhAEdnd/5ANOcxwgr3Pws1L8Z/vs4JG+1Oo0oQ47cXJIeHkzu9u3UfHUSlf7zb6sjCSEsYvYKf6Tz2A+AtlrrdVrrMw+Jtdb5Wuu1GFf4KzCK/mOlDStMsvkY8/PDqhuL8mSlWZ1IlAFHfj7JQ4aS8+OP1Jj4CpU6dbI6khDCQmYL/k2ABkZorS+4RZvzu+HOtm1N9iVcIbgK9FhkbKO7/H6wF1qdSLiRLiggZfgIY4vbceOofPvtVkcSQljMbMGvCpzUWh++VEOt9SGM7XSrmuxLuErMFXDHVDi4GdY+Z3Ua4Sa6qIiUUU+StWED0c8/R3j3blZHEkJ4ALMF/xQQVpwNc5xtKjmPEVa74h645mH4YSb8vMzqNMLFtMPB4WeeIfPzz6n25JNU6d3b6khCCA9htuD/hPFc/tFitB3mbLvNZF/C1f49Fuq0gc8ehcM/W51GuIjWmtSXxpDxyadUHfYokf3utzqSEMKDmC34szCm3Y1VSo1TSlU+t4FSKkYp9TowBuMZ/izzMYVL+fjB3fMhOBKW9oHsY1YnEqWktSZt0qucXLaMyEGDiHr4YasjCSE8jKmCr7X+CFjoPH40kKqU+p9S6kOl1GqlVCJwAOPq3gYs0Fp/7KrQwgVCqxrb6WYdgRX9wC67F5dn6dNncHzePCL69KHqiOFWxxFCeKDSrLR3P/A0kAkEAK2BrsCtQDPA3/ndU8CAUqUU7lGzFdz2OhzYBOtfsjqNMOnY3HmkT5tG5W7diH56NEopqyMJITyQ2e1x0Vpr4BWl1FvALcBV/D0S/yjGc/61WuucUqcU7nNlHzi0Hb6bCjXioHl3qxOJEjix7APSJk0i7NaOxIwdg7LJflhCiPMzXfBP01pnAyudL1Ee/WcCpO6CT4ZAVCOo3tzqRKIYMlavJvXFFwm96SZqTpyI8vGxOpIQwoPJ5YAAX3+4ZwEEVoalvSHnuNWJxCVkbtzIoSefIrhVK2pOeRPl7291JCGEh5OCLwxh0XDPQsg8DB8OAIfd6kTiArK3bCFl2HACGzUi9u2Z2AIDrY4khCgHLnlLXyn1lfPHg1rrfud8VhJaa93exHGirNS6Gjq9ZszPX/8S3DLG6kTiHLmJu0h+eDB+sbHUmjMbn9BQqyMJIcqJ4jzDb+t8/+U8n5WENnGMKGut7oPDO2DzFIhpKYP4PEj+/v0kDRqET3g4tee+i29EhNWRhBDlSHEKfj/ne8Z5PhPeqONEOLLbOYivIVRvYXWiCq/w0CH+6j8AfHyoPfdd/KKjrY4khChnlDG7zjvFx8frrVtl/3dTMlNhVlvw8YdBG43d9oQlio4f5+C9fShKT6fOwgUENm5sdSQhhIdSSm3TWsef7zsZtCfOL6z634P4ZCU+y9izskga+ACFhw9T6+2ZUuyFEKaZKvhKqa+UUstL0H6JUmq9mb6EhWpdDZ0nw/6NsP5Fq9NUOI78fJIfGULevn3ETnmT4FatrI4khCjHzC680xZILUH7a4HaJvsSVrqqr7Gj3ndvQfWWcMXdVieqELTdzqHHnyDnhx+o8eokQm+6yepIQohyrqxu6fsgo/TLr/9MgNr/gk+HyHa6ZUBrTeqLL5G5bh3RT4+m8u23Wx1JCOEF3F7wlVIBQDXglLv7Em7i6w/3vOfcTvdeyE63OpFXOzplCieXLyfywQep0rev1XGEEF6iWLf0lVK1gbrnfOyvlLoBuNDWXAoIB3ph7Jz3ncmMwhOEVoMe78PcjrD8fkj4GHz8rE7ldY4vWMixt98h/O67qTp8mNVxhBBepLjP8PsBz5/zWQSwsRjHnv4HwZvF7Et4qppXwR1T4eMHYe2zcOtEqxN5lYzVqzkyYQKhHdpT/cUXZJtbIYRLFbfgnwT+Ouv3OoADSL7IMQ6M2/i7gXe11htMJRSepWVPOLwT/jfdWJDnyj5WJ/IK2d99x6GnRhub4UyeLDvfCSFcrlgFX2s9BZhy+nellAM4qrWu565gwoPdMgaO7IJVI6BqY4g97xoPophyd+8mechQAurVI3bGdGwBAVZHEkJ4IbOD9l4CJrsyiChHfHzh7vkQFgPL+hir8glTCg4eJGnQg/iEh1Nr9mx8KlWyOpIQwkuZmoevtX7J1UFEORNcBXouhndvgWUJcP8q8JUr05IoOnaMvx4YBHY7tebMwS+6mtWRhPAqhYWFJCcnk5eXZ3UUlwsMDCQ2NhY/v+IPnja78I4QUL05dH0bPugLq0fCHW+BDDQrFkd2NkkPPkTR0aPUmT+PgMvk6ZgQrpacnExYWBh169b1qkGwWmuOHTtGcnIy9eoV/++OUhV8pVQQcBdwPVADCOHC0/S01rp9afoTHqjpnXDjE/D1q1D9CrhmkNWJPJ4uLCR52HDy9u4ldvo0glq2tDqSEF4pLy/P64o9gFKKyMhIjh49WqLjTBd8pdTNwGKgKkaRP72S3tn/y579WalX2lNKhQNzgObO8/XXWn9f2vOKUmr7NKTugs+fgmpNoN4NVifyWFprDj/7HNnffkvM+HGEtW1rdSQhvJq3FfvTzPx3md0853LgE4wV9NYDIzCK+ilgIPAMsMH52TFgKNDfTF/nmAJ8rrVuDLQE9rrgnKK0bDboNgsiLzdu7584aHUij3X0jTfJ+OQTqg57lPDu3a2OI4SoQMyO0n8C4/b9+1rrfzun7QHkaq3naq0nOG/fdwQCMRbuWVqaoEqpysCNwLsAWusCrfXJ0pxTuFBgJei1BLQdlvaGgmyrE3mc44sXc2zWLMJ79CDyoYesjiOEqGDMFvybMW6pj7tYI631WmA4cBXwuMm+TqsHHAXmKaW2K6XmKKVCzm2klBqklNqqlNpa0ucbopQi68NdcyFtD6x8GLTsl3Ra5pdfcmTsOELbtaP6c8967W1GIYTnMlvwawIFWut9Z33mwLiaP9dioAi4x2Rfp/li/MNhptb6SiAbeOrcRlrrWVrreK11fNWqVUvZpSixyzsYC/Ps+QS+fs3qNB4h56ftpIx8nMArWlDz9ckoX5kcI0RFNnDgQFatWsU777xDTEwMcXFxZ16JiYlu69dswc8Hss75LBOorJTyP/tDrXUeRnEu7byjZCBZa/2D8/cVGP8AEJ7muiFwRQ/YMA5+WW11Gkvl7z9A8sMP41s9mlozZ2ILCrI6khDCYtu3bz9T3MeNG8eOHTvOvFq0aOG2fs0W/GSM4n72pcofzvf/t86qUqo6UJkLT9crFq11KpCklGrk/Kg9sKc05xRuohTcPgVqXAkfDYK0ijm2sig9naRBg8DHh9qzZ+NbpYrVkYQQFti3bx9t2rShRYsWjB8/ntTUVGJjY9m5cydxcXFllsPsvcU9QFOMkfLbnJ+tx7jifl4p1UVrnee82j89oG97qZIahgKLnOfdjzEYUHgivyBjJb5ZbWFJT3hgg7E6XwXhyMkh6aGHKUpPp87CBfjXrm11JCEqtJc+282eQ6dces6mNSrxwu3NLtomPz+frl27Mm/ePFq3bs3gwYNp3LgxALt376Zfv37YbMa19+DBgxk0yH1rmZgt+GuAu4E7+bvgTwUeBm7BuBLfBzQAIjEG+JV67X2t9Q7OuYMgPFilGtBjEczvBMvvgz4fgU/xl4Esr3RRESmPjSRvzx5ip00jyI236IQQnm3lypXEx8fTunVrAJo1a0ZgYCBJSUlUrVqVnTt3llkWswV/Bca0vPTTH2itU5RStwPvYwzqu875VQ4wWmv9SWmCinKq1tXG7f2VD8MXz0CnSVYnciutNanjx5O1cSPVX3iesJvbWR1JCAGXvBJ3l8TERFq1anXm923bttG2bVsSExNp0qRJmWYxu3lOFjD9PJ9vUkrVwyj2sUAGsFlrnVGqlKJ8i+ttrMT3v+kQ3Qxa3Wd1Irc5NmcOJ5csJfKBgUT06mV1HCGExSIjI9m1axdgFPslS5YwbNgw1qxZcz4u+R0AACAASURBVObWflkxVfCVUqf38MzWWtvP/k5rXQR8U9pgwsvcMgaO7jU22YlqCHWuu/Qx5UzG6tUcnfw6lTp3puqIEVbHEUJ4gISEBDp16kRcXByNGjUiPDycpk2bMmnSJDZt2sSaNWsAY6ncb775htDQULdlMXtL/yTGvPt6QJLr4giv5eNrLMozuz0s6wODNkC49wxky/nxRw4/NZrg+HhiJryMspmdACOE8CZRUVFs2bLlH58vWrSozLOY/VspCziltZZiL4ovKAJ6LQV7ISzpDfnnLuVQPuXv30/SkKH4xcYSO+0tbP7+lz5ICCHKmNmCfwAIPmcevhCXVrWhc/nd3cZAPofD6kSlYsy1fxDl60ut2bPwCQ+3OpIQQpyX2YL/AeAHdHFhFlFRNOgAt4yFvZ/CpolWpzHNkZtL0uBHKEpPp9bbM/GPjbU6khBCXJDZgv8qsBV4RynV3oV5REVx3SMQdy9segV2f2x1mhLTDgeHRj1JXmIiNV97VebaCyE8ntlb8k8BXwFNgLVKqZ3A9xi72dkvdJDWeozJ/oS3UQpuewPSf4OPH4Yql0FMS6tTFVvaq6+RuW4d0aOfIqxDB6vjCCHEJZkt+C9irJ53en38lsAVF2mvnO2l4Iu/+QZAz0Uwqx0s6WUsvxsWbXWqSzqxZAnH580j4t57iejb1+o4QghRLGYL/gKMAi5E6YRWg16LYW5HWHYv3LcK/M63y7JnyNq0idSx4wht25bop0fLvvZCiHLD7Ep797s4h6jIYlpC17fhg76wajh0mWnc8vcweb/8QsqIxwho3Iiak19D+fhYHUkIIYpNVgcRnqHpndDuGfh5CWyecun2ZazwyBGSHnwIW6VK1Jr5NraQEKsjCSHKqYEDB7Jq1SoOHz5Mz549iY+Pp2HDhrRr5969N0wVfKXUXKXU6yVoP0kp9a6ZvkQFcuMT0KwbfPki/PJfq9Oc4cjOJunhh3FkZlLr7Zn4RVezOpIQohzbvn07cXFxJCQk0LVrV7Zu3cq+ffuYOnWqW/s1e4V/P9CzBO3vdh4jxIUpBV1mQI04+OgBOLLb6kRou52Ux58g/5dfqfnG6wSW8WYXQojyb9++fbRp04YWLVowfvx4UlNTiYmJYePGjdx0001n2rVw8/Teslop7/QofSEuzi8Iei6B2e1gcU944CsIrWpZnCMTJ5K1YQPRzz9H6Fl/MIUQ5cyapyA10bXnrN4Cbn3lok3y8/Pp2rUr8+bNo3Xr1gwePJjGjRvj4+NDhw4daNmyJV26dKFv375cf/31rs13Drc/w1dK2YBqQLa7+xJeolKMMV0vO83YaKco35IYxxct4sSChVS5ry9Veve2JIMQonxbuXIl8fHxtG7dGoBmzZrRsqWx5siaNWv48MMPqVy5Mh07dmTlypVuzVKsK3zndrjnLhLuo5Sqxd9z8f9xmPOYvkAg8LPZkKICqtnKuL2/oj+sGgF3Ti/TkftZ33zDkfEvE9q2LdVGjSqzfoUQbnKJK3F3SUxMpFWrVmd+37ZtG23btgWMLXHbtGlDmzZtOHHiBDt37qRLF/etWF/cK/wRGBvmnH4BRAF/nvP52a/9wE/AcGf72S5JLCqO5t3hpqdgxyL4zr2DWc6Wt28fKcNHENBIpt8JIUonMjKSXbt2AUaxX7JkCS1btuSLL76goKAAgLS0NL799ltuueUWt2Yp7jN8xf+/ktdc+Mr+7DangN3AHK31/BKnE+KmJyH9V1j3AkQ2gMad3Npd0dGjJD30ELbgYGrNnCHT74QQpZKQkECnTp2Ii4ujUaNGhIeH07RpU6ZNm8bgwYMJDQ0lICCAsWPHct1117k1S7EKvtb6RYzldAFQSjmAVK11DffEEsLJZjMW4jlxED4cCAPWQvXmbunKkZdH0pAh2E+cpM7ChfhVr+6WfoQQFUdUVBRbtmz5x+ezZ5f9TW+zg/YWYGyRK4T7+QVBz8UQWBmW9ISsNJd3oR0ODo0eTd7ORGpMmkhQ82Yu70MIIaxkquBrre/XWg+/dEshXKRSDPRaAjnHYGlvKMx16enTp00nc83nVBv5GJXc/BxNCCGsIEvrivKjRhx0fQeSf4RPhoB2zdIOGZ99RvqMGVTu3o0qAwa45JxCCOFpSrXwjlKqCdAdaA5EAH4Xaa611u1L058QNL0D2r8A61+CqIbQ9slSnS7np584/PQzBF99NTEvvCC73wkhvJbpgu9cS/9R/jmC/0JkpT3hGm1GQPo+2PgyRNaHFneZOk1BcgrJQ4biV6MGNadOQfn7uzioEEJ4DlMFXyn1CH/Pr08EPgFSgDwX5RLiwpSC26fAiT9h5WAIrwO1ri7RKexZWSQ//BC6qIjYmTPxjYhwT1YhhPAQZq/wH8C4Yn9LBu8JS/gGQI/3YU57WNrLWHM/vHaxDtV2OykjR5K//wC158wm4LJ6bg4rhBDWMztor6Hz/XlXBRGixEKioPcHUFQAi3tA3qliHZY2aRLZm76m+nPPEuLmhS6EEMJTmC342UCG1rp4f8MK4S5VG8E978HRX4119+1FF21+YtkHHH9vAREJCUT0LMkOz0II4RoDBw5k1apVvPPOO8TExBAXF3fmlZiYyOHDh+nZsyfx8fE0bNiQdu3auaRfswX/B6CSUsq6fUuFOK1+O+g8GX5fB188fcFm2f/7gdSxYwm54Qain5QNcYQQ1ti+ffuZ4j5u3Dh27Nhx5tWiRQsSEhLo2rUrW7duZd++fUyd6pq9RMwW/AkYz/CfcUkKIUorvh9cNwS2vAM/zPrH1wUHD5I8bBj+depQ8/XJKN9SzUgVQohi27dvH23atKFFixaMHz+e1NRUYmNj2blzJ3Fxcf+vrd1uZ+PGjdx0001nPmvRooVLcpj6W09rvVkpNRB4WykVCLyitf7TJYmEMOuWMXB8P3z+JETUhYb/BsB+6hRJDw9GKUWtmTPwCQuzNqcQosxN3DKRX47/4tJzNq7SmCdbX3wtkPz8fLp27cq8efNo3bo1gwcPpnHjxgDs3r2bfv36YbMZ196DBw9m0KBBdOjQgZYtW9KlSxf69u3L9ddf75K8Zqfl7Xf+aMcYsf+AUuo4kHmRw7TWur6Z/oQoFpsPdJsN826FFf2g/xfoqMakjHiMgqQk6sx9F//axRvJL4QQrrBy5Uri4+Np3bo1AM2aNSMwMJCkpCSqVq3Kzp07/3HMmjVr2Lx5M59++ikdO3Zk4cKFdOnSpdRZzN7XrHuezyKdrwuRhXeE+wWEQu9lMLs9LO7BkYwuZG/eTMy4sQRfXbK5+kII73GpK3F3SUxMpFWrVmd+37ZtG23btiUxMZEmTZqc9xilFG3atKFNmzacOHGCnTt3WlrwXTNkUAh3qFQDei/lxFN3cuKHj6iS0Jvwu8ytxieEEKURGRnJrl27AKPYL1myhGHDhrFmzZozt/bP9sUXX9CuXTv8/f1JS0vj22+/Ze7cuS7JYvYZ/iaX9C6Em2T/lU/q1hBCYnKpVm8vOOzGLX8hhChDCQkJdOrUibi4OBo1akR4eDhNmzZl0qRJbNq0iTVr1gDGVf0333zDihUrGDx4MKGhoQQEBDB27Fiuc9F6ITJUWXidgr/+IuXRR/GvW4+aI29Bff08rHse/jPe6mhCiAomKiqKLVu2/OPzRYsWnbf97Nmz3Zal1AVfKeULtAJqAcFa6wWlTiWESfasLJIeHgxArRnT8alTB/JS4Ptpxsj91g9YG1AIISxS2u1xnwSewNga97QFZ30fDnwH+AM3aq0PlaY/IS7m9Br5BQcPUnvOHPzr1DG+6DgBTh6ENaOMjXac0/WEEKIiMbvwDkqpRcDLGMX+APCPNU211ieBTUA9QNYxFW6VNvl1Y438Z54m5Npr/v7C5gPd34Xo5sZ0vdRE60IKIYRFTBV8pVRPoBeQCvxLa305cPwCzRcBCuhgKqEQxXDy45UcnzuXiN69iOjV658NTk/XC6wMi+6BjJSyDymEEBYye4U/AGNe/XCt9Q+XaLsVcADNTfYlxEXlbN9O6vPPE3zNNUSPHn3hhpVqGLvr5Wcau+vlX2ydKCGE8C5mC/6VGAX/00s11FrnARmAbLQjXK7w8GGShz6Kb0wMNd98A+Xnd/EDqjeHe+ZD2h5Yfv8ld9cTQghvYbbghwKZWuv8Yrb3x1iGVwiXceTmkvzIEHRuLrVmTMc3IuLSBwFc3gFuex1+/xL+OxK0LAIphPB+ZkfpHwVqKKXCtNYXvS+qlGoAhAD7TPYlxD9orTn09NPk7d1L7MwZBFx+eclO0Op+OHEQvn3dmK7XZoQ7YgohhMcwe4W/2fl+dzHaPoFx+3+Dyb6E+Idj77xD5prPqfb4SMLatjV3kpufg+bd4csXIXGFK+MJIcQFDRw4kFWrVvHOO+9QvXp14uLiuOyyy5g/f75b+zVb8N/CGHk/Til13sF4SqkApdR4YCBGwZ9msi8h/p/M9es5+uYUKt1xO1X69zd/IpsNusyEOtfDyofhz82XPkYIIUpp+/btxMXFkZiYyIsvvsiOHTtYsWIFI0eOdGu/pgq+1noz8CpQHfhBKfUJEAaglHpdKbUCSAGech7yvNZ6twvyigoub98+Dj0xisAWLYgZMwalVOlO6BsAPd43busv7QVHf3VJTiGEOG3fvn20adOGFi1aMH78eFJTU4mNjWXnzp1nNtCJjY3FbnfvUDfTK+1prZ9USh0CxgK3n/XVMIyrf4BsYLTWWq7uRakVnThB8uBHsIWEEDvtLWyBga45cXAVuHc5zOkAi+6CAV9CWLRrzi2E8AipL79M/t5fXHrOgCaNqf700xdtk5+fT9euXZk3bx6tW7dm8ODBZ4r86S1ytdZMnTqV2267zaX5zlWqpXW11lOUUvOB7sC/gBiMuwZHgO+B5VrrCy3IY4pSygdjbn+K1tq9/+sIj6ELC0kZPoKitDTqLFyAX7SLC3JEXWOO/vzOsPgeuH+1sViPEEKUwsqVK4mPj6d169YANGvWjMDAQJKSksjKyuI///kPfn5+tG7dmunTp7s1S6k3z9FaZwBzna+yMAzYC1Qqo/6EBzjyykRyfviBmFcmENSypXs6qXkV3DXPuLW//H7otRR8ZENJIbzBpa7E3SUxMZFWrVqd+X3btm20bduWxMRE2rdvz+eff15mWUyvpW8FpVQs0BmYY3UWUXZOLF/OiUWLqHL//YR36eLezhp1hM6T4fd1sHqEzNEXQpRKZGQku3btAoxiv2TJElq2bMnOnTtp6a6Llwswu5Z+FaVUX6XU7cVoe4ezbbiZvs7xJjAKY6neC/U3SCm1VSm19ejRoy7oUlgpZ9s2UseMJaRNG6o97t4RrGfE94cbHoefFsDXr5ZNn0IIr5SQkMCOHTuIi4tj0qRJhIeH07RpUxITE7niiivKNIvZ+5V9gcnAGOCzS7S9ERiBcSve9OA9pdRtQJrWeptSqu2F2mmtZwGzAOLj4+XyrBwrPHyY5EeH4V+jBjUnv4byLcPb6zc/C6dSYMN4qFQTrry37PoWQniNqKgotmzZ8o/PFy1aVOZZzN7S7+p8X1qMtnMxRu13M9nXadcDdyil/nT2e7NS6v1SnlN4qDPL5ublETtjOj6VK5dtAKXg9qlwWVv4dCj89mXZ9i+EEC5mtuDXB/K01pectKy13gPkASVc+/Qf5xmttY7VWtcFegJfaa37lOacwjNprTn87HPk7d1LjddeJaB+fWuC+PrDPQshuil80BcObbcmhxBCuIDZgh8F5JagfQ5QzWRfooI5Pncup1avpurw4YS1a2dtmMBKcO8KCI6ERXfD8QPW5hFClIj20oG3Zv67zBb8E0C4UuqSU+OUUpWBcIwtcl1Ca71R5uB7p6xvviHttcmE3dqRyEEPWB3HEFYd+nwIjiJ4vztkp1udSAhRDIGBgRw7dszrir7WmmPHjhFYwsXHzI6C2gbcCiQAl1opIAHjHxY7TPYlKoj8AwdIeWwkAY0aUWP8+NIvm+tKVRtCr2Ww4A5jYZ77PgP/EKtTCSEuIjY2luTkZLxxxlZgYCCxsbElOsZswV8IdAImKqX2aa3Xna+RUurfwCsYm+csNNmXqADsWVkkPzIE5etL7LRp2IKDrY70T7WvgbvmwrI+xsI8PReDj5/VqYQQF+Dn50e9evWsjuExzG6eswxYDwQDa5RSq5RSjyilbne+hiil/gv819lmk9ZaRtSL89IOB4eeGEXBwYPUfPNN/GNrWh3pwhp3hs6vw29r4dNHZWEeIUS5UZqJzd2BxRhX+p0wbvGf7fT92DVA71L0I7xc+rRpZG3YQPSzzxJyTWur41xafD/IOgIbJxjP9zu8YHUiIYS4pNLslncKuE0pdSvGQjzXAqd3NDkC/A9YoLVeU+qUwmudWruW9Bkzqdy9GxH3lqN/F970JGQehm9fN4r+NQ9anUgIIS7KFZvnrMG4iheiRPJ+3cehp0YT1LIl1V94wbMG6V2KUtBpMmQdhTVPQkgUNO9udSohhLigcrV5jvAe9pMnSR4yBJ+QEGpOnYrN39/qSCXn4wt3vQu1r4OPHoQ/vrI6kRBCXJAUfFHmdFERKY+NpCg1ldi3puIXXY7XZPILgl5LIKohLEuAlJ+sTiSEEOdVqlv6zh3wbgOaAxHAxeYoaa31gNL0J7xD2htvkP3dd8SMG0tQXJzVcUovKNxYmGfuv43V+AashUiLlgMWQogLMF3wlVKPAhOA00v9XOoBrAak4FdwGatWc/zduUT07k34XXdZHcd1KsVAn4+Nor+wC/Rfa3wmhBAewlTBV0r1xNibHuAo8AWQgrFJjhDnlbdnD4effZbg+HiiRz9ldRzXi7rcWHf/vdvh/W5w/2oIrmJ1KiGEAMxf4Q9zvi8H7tNaS6EXF1V04gTJQ4biEx5OzSlvovy8dIW6mldBz0XGrf3FPaDvSlmCVwjhEcwO2muOcYt+iBR7cSm6qIiU4SMoSk8n9q238I2MtDqSe13WFrq/CylbjW11iwqsTiSEEKYLfhGQobX2vh0JhMulvfoqOT/8QPUxLxHUornVccpG0zvgtjfg9y9h5cPgcFidSAhRwZm9pb8DaKOUquRccU+I88r45BOOv7eAiIQEwrt0sTpO2Wp1P+Qch/UvGSP5O71mLNgjhBAWMHuF/zrgAzziwizCy+Tu3s3h518g+OqriR71hNVxrNFmBPzrUfhxDmwYb3UaIUQFZuoKX2v9mVLqeeAlpZQGpmitc10bTZRnRcePkzx0KD5VqlDzzTe8d5DepSgFt4yBvJPw9asQGA7/GmJ1KiFEBWR2Wt7pNUSzgfHAc0qpPUDmRQ7TWuv2ZvoT5YsuKiJlxGPY049RZ/Fi7x+kdylKwW1vQt4pWPsMBFaGqxKsTiWEqGDMPsNve87vQUCrSxwjG4dXEGmvvkbODz8Q88oEgpo3szqOZ7D5QLdZkH8KPnsUAsKgWQUb0yCEsJTZgt/PpSmE18j47DOOv/dexRykdym+AdDjfVjYDT4cCP6h0KCD1amEEBWE0tp7L7zj4+P11q1brY5RYeTt2cOfve8lqEULas99t+I+t7+U3JPGanzpvxlr8Ne93upEQggvoZTaprWOP993sluecIn/t5JeRR6kVxxB4ZDwMYTXMlbjkx32hBBlwGUFXykVpJSq5XwFueq8wvMZ290+5lxJb6oM0iuOkChIWAnBEca6+0f2WJ1ICOHlSlXwlVJVlFIvnjVC/0/nK1MptUcp9YJSKqL0MYUnS3vjDXK+/x/VX3iBoBYtrI5TflSuCX0/Ad9AWHAnpP9udSIhhBczXfCVUq2BXcBzQGPnuZTzZXN+9jywy9lWeKFTa9Y4t7vtRXj3blbHKX+qXGYUfe2ABXfAiYNWJxJCeClTBV8pFQ2sAaoDJ4EJwC1AE+frFuAV53cxwGrnMcKL5P26j0NPP0PQVVcR/ZQXbndbVqo2Mp7pF2QZRf/UYasTCSG8kNkr/FFABLATaKK1fkZrvV5r/avztV5r/TTQFEgEqgAVdG1V72TPyCB5yBB8QkONQXr+/lZHKt9iroA+H0F2ulH0s2RfKiGEa5kt+J0xFtLpr7VOu1AjrfURoD/Gbf7bTPYlPIx2OEh54gkKU1OpOXUKftWqWR3JO8TGQ+9lcDLJeKafc9zqREIIL2K24NcGMrXWl5xPpLXehjGgr7bJvoSHSZ82jeyvv6H6M08TfOWVVsfxLnXbQK8lcOx3WNjFmLMvhBAuYLbgFwD+Sl16r0+llA3wcx4jyrnM9etJnzGTyt26Ed6jh9VxvFP9dsaKfEf2wPvdjTX4hRCilMwW/F+AAKBrMdp2BQKBX032JTxE/v4DHBr1JIHNm1P9hecpxr/3hFkN/w13z4fDO2DxPZCfZXUiIUQ5Z7bgf4DxXH6WUuqCi4Erpe4AZmE8719qsi/hAexZ2SQPHYry9yd26hRsAQFWR/J+TW6D7nMg6QdjRb6CbKsTCSHKMbOb50wD+gBxwBdKqa3ABiAF42q+NnAT0AzjHwbbgemlTissobXm8OjRFBw4QO257+JXo4bVkSqOZl3BYYePHjCKfu8PwD/Y6lRCiHLIVMHXWhcopf4NLAT+A1wNnLtY/+n7vZ8DfbXW8gy/nDo2Zw6Z69ZRbdQoQq691uo4FU+Lu4yi//GDsLQX9FoKfrJ6tRCiZMxe4aO1TgduVUq1Ae4CrgKqOr8+CvwErNBaf1vqlMIyWZs3c/SNNwm7tSNV+t1vdZyKq2UP0HZYORiW9oaeS8Av0OpUQohyxHTBP81Z0KWoe6GC5BQOPTaSgPr1qTFunAzSs1pcb+NK/9OhzqK/WIq+EKLYil3wlVJtgRsx5t+/UcxjHgNCga/kSr98ceTlkfLoo2iHg9i3pmILCbE6kgC4KsF4/3SocXu/52K5vS+EKJZijdJXSgUC7wMvAMklOH8S8CLwnlJKNkgvJ7TWpL74Enl79lBj4kT869a1OpI421UJcOc0+GMDLOkFhblWJxJClAPFnZbXHagBbNRaLy/uyZ1tNwF1Kd6cfeEBTi5dSsbKlUQNHkzYze2sjiPO58o+cOd02L9Rir4QoliKW/C7YMyln2aij7cwRux3N3GsKGM527eT+vIEQm68gaghj1gdR1zMlfdClxlG0V98j8zTF0JcVHELfivn+5cm+jh9zLnT9oSHKUpPJ2XYcPyqV6fmpEkom9l1mUSZiesNXd+BP7+F9++C/EyrEwkhPFRx/0avhjFYr8R/m2itT2FsnhNd0mNF2dGFhaSMeAz7qVPEvjUVn/BwqyOJ4mrZ4+8V+RZ2g7wMqxMJITxQcQu+KkHbCx0vc7o8WNprk8n58UdixrxEYOPGVscRJdW8u7H2/qGfYEEXyD1hdSIhhIcpbhFPB0KUUhEl7cB5TKjzHMIDnfrvfzn+3ntE9OlD5TvusDqOMKvpHc5d9nbBe7dDtvyRE0L8rbgF/2fne0cTfXRyvu80caxws/zffuPQs88RdOWVRI96wuo4orQa3Qq9lkD6bzC/M2SmWp1ICOEhilvw12Dckn9WKVXsbdKcbZ/BGOG/uuTxhDvZMzNJHjIUW0gwNd98E+Xvb3Uk4QqXd4A+H0JGMsy7FU4mWZ1ICOEBilvw5wOpQGNguVIq7FIHKKVCgeXOY9KA90xmFG6gHQ4OPTWagpQUYt94A7/oalZHEq5Utw0krIScY0bRP/aH1YmEEBYrVsHXWucCD2JcqXcGdiulRiqlGp7bVinVUCn1OLDb2dYBPOg8h/AQx2bPIWv9eqJHPUFwvMyY9Eq1rob7PoPCHJjbEY7stjqREMJCxR55r7X+DBgIFACxwCRgr1IqRymV4nzlAHuBiUAtZ9tBWutPXR9dmJW1eTNHp0yhUufORCQkWB1HuFNMS+i3Bmy+MK8TJP1odSIhhEVKNNVOaz0fuIa/n8crIBCIcb4C+Xv63WrgWq31XJckFS5RmJLCoZGPE1C/PjFjx8gOeBVB1UbQ/3MIioAFdxor8wkhKpwSb4+rtd4J3K6UqgG0BZoAkc6vj2Fc4W/SWqe4KiSAUqoWsABjAR8NzNJaT3FlH97OkZ9P8rDh6KIiak6dgi042OpIoqxE1DGK/sKusOhuuGsuNLnd6lRCiDJU4oJ/mtb6ELDYhVkupQgYqbX+yTlocJtSap3Wek8ZZijXjowbT96uXcROn0ZAvXpWxxFlLaw63L/aKPgf9IXbp8BVfa1OJYQoI+VmsXSt9WGt9U/OnzMx7iTUtDZV+XFyxQpOLl9O5AMPENa+vdVxhFWCq8B9n8Jl7eDTofDN66C11amEEGWg3BT8syml6gJXAj9Ym6R8yN21m9QxYwn513VUHT7M6jjCav4h0GspNL8L1r8Ea58Fh8PqVEIINzN9S98qzvn9HwLDnRvznPv9IGAQQO3atcs4necpOnGClGHD8ImMpMZrr6F8fKyOJDyBrz90mw3BkfD9NMhKgzunG58LIbxSuSr4Sik/jGK/SGv90fnaaK1nAbMA4uPjK/S9Sm23c2jUkxSlpVFn0fv4VqlidSThSWw2uHUihFaFr8ZB9lHosRACLrmulhCiHCo3t/SVMX/sXWCv1vp1q/OUB+kzZpL9zTdEP/MMQVdcYXUc4YmUghufMK7uD3ztXH//iNWphBBuUG4KPnA9kADcrJTa4Xx1utRBFVXWpk2kT5/O/7V37/GRVGXCx3+nLn3v3DOZJJO5ZBhAUBAYWERBFEHlVVBQwdcLuuJ1EdQVFtQXEVwVdpWLIooru7iKIqgsrrgKAiugCIwKAgqTuZLM5Da59r276rx/VHUnmUlmkiGT7k6e7+dzPqfq1OnOU+nLU1Vddar2LW+h7px3lDscUemOerf3u/7gRvjuqTDYVe6IhBDzTOlFfIbu+vXr9RNPPFHuMBZcrrubLWe/DbutjdU/vA0jFCp3SKJadG+A297unbn/zh/CyuPLuvA48QAAIABJREFUHZEQYg6UUhu01tOOl15Ne/hiFtxMhu4LLwStWXHD9ZLsxdysOAY+cK83Kt+tZ8DT054qI4SoQpLwF5neq64i++xfabv6KwQ6OsodjqhGjWvh/Pug7Si48/3w8HVyrb4Qi4Ak/EVk+I47GP3JT2n86EeIv+Y15Q5HVLNIA7z3v+Dwt8J9n4f//iQ4+XJHJYR4Earqsjwxs/TTz9B31ReJvvKVNF9wQbnDEYuBHYKzb4G6VfDIdTC8Fd7+HxCuK3dkQoj9IHv4i4AzMkLPhRdiNjXS9q//IoPriPljGHDqF7zL9rY+7J3BP7S53FEJIfaDJPwqp12XnksuoTAwwIrrr8eqry93SGIxOurd8N67vMF5vnMKbPtduSMSQsyRJPwqN3jTTSR/+xAtn/0M4Ze9rNzhiMVs9avg/N94v+/fegZsuLXcEQkh5kASfhVLPPQQg9+4kdozz6DunHPKHY5YCopn8K85CX5+IdxzCTiFckclhJgFSfhVKt/Tw45PX0xw3TqWX3EF3sjDQiyAcD383x/DKy6Ax74N3z8LUkPljkoIsQ+S8KuQm8vRfdEn0I7jDa4TDpc7JLHUmBa8/p/hzG/C9t/DzSdD71/KHZUQYi8k4Vehvi99iczTT9P2lS8TWL263OGIpeyod8H77gEnB/92KvzlznJHJISYgST8KjNy112M/Oh2Gs//APHXva7c4QgBHcfCh/7XG5nvJx+A//mM/K4vRAWShF9FMs89R+8VXyBy7LE0f+IT5Q5HiAnxFjjvbjjuw/DojfC9M2G8t9xRCSEmkYRfJZzxcbovvBAzHqf9a19FWTJIoqgwpg2nXwNvvRl6NsC3ToQtD5U7KiGETxJ+FdBas+Oyy8h399B+3bVYzc3lDkmImR15DnzwfgjVwvfOgIe+Bq5b7qiEWPIk4VeBoVtuIXHfb1h28aeJHHNMucMRYt9aDoMPPQCHvQV+8wX44bmQ3FXuqIRY0iThV7jkY4/R/7Vrib/+9TScd165wxFi9oJxeNstcPq/wuYH4Fuv9MbjF0KUhST8Cpbv76fnU/9IoKOD1n/+ogyuI6qPUnDcB70heQNRuPXN8MCXwXXKHZkQS44k/Aql83l6PvUp3GSS9huux4zFyh2SEPuv9Qjv0r0jzoH//Qr8x5tgZHu5oxJiSZGEX6H6r72O9BMbaL3yC4QOPrjc4Qjx4gVj8NZvwVu/7Y3Kd9Mr4akfg9bljkyIJUESfgUau/dehm65hbp3nkvtm99c7nCEmF9HngsffRiWHQY//aA3WE96uNxRCbHoScKvMLmtW9l52WcIvexltFx2WbnDEeLAqF8N7/sFvPZz8Ox/eXv7XfeVOyohFjVJ+BXETafpvugTKNNkxXXXYgQC5Q5JiAPHtOCki+ED90IgBt8/G+6+ELLj5Y5MiEVJhmubpUe6BukZSfOO9R0H5Pm11vReeRXZ55+n4+ZvY7e3H5C/I0TFaT8aPvxbePBL8Luvw6YH4MyvQ+fJ5Y5MLBJaaxxXU3C92tEax/HmXX+ZM2na1Uya1mhNadrV3vMV+xSX68nz/t/02vxp8PtNmkazvCbE+tUNC/J/kIQ/S3f9qYcHnhs4YAl/5M47Gf3Zz2j62MeInXjiAfkbQlQsOwSnXgmHvgnu+qg3Fv9R74HTroJwfbmjEy9C3nFJZR2SuQKpnEM655DOO6RyBTJ5h0zeJZ332rMFl2zBa8sWvPlcwfVrh1zBJe9ocgWXnOOS90vB0eT8uuB6fQqOW0rwBbdyTww97bAWSfiV5pDlce7Y0M2uRJbGWHBenzvz7LP0XfVFoiecQNM/fGxen1uIqtJxHHzkYXjwK97e/vO/gtP/BQ4707umXyworTWJbIGRVJ7hVI7RdH5KGUsXGM/kGc8UGMvkSWYLjGcKJHMFEpkCyZyXpOcqYBoELYOA5dVB2yRgevMBy8A2FXHbImAa2KaBZSoCfm2ZBrbh1ZahsEyFaXhthqGwDIU5qTYMhakmatNQKAXmpHZDKQyFVxsT0wowjKm18vuq4nKlSm9dpSbmFV4dCy5cGpaEP0uHLI8D8FzfOCfMkPBd7TKaHWU4O8xIZoREPrFHH1OZGMrAVCamYWImMgQu+CzUxUl99sNsHt+KoQwsZXnLlYllWFiGVXpMcZmh5BQMsQjZYTj1C/DSs+Duj8Md58Ehp8Mbr4a6leWOruqlcgX6x7IMJLIMjmcZTGQZTOQYTGQZTuXYlcgxnMoxlMwxksrvde/YNhXxkE08ZHklaLOiPkI8ZBENmkSDFtGARTRoEQmYRAImYdskErAIB0xCtkHYNgn77UHLJGgZGIZs3B0IkvBn6ZAWL+E/3zvOCWubpiwbzY5yzePX8IvNv8DRsx9BTGnNp3/iclSv5vPvMtn48AfmFFNxw8EyrNJGQHHDoLiRMLm9OF+anjRvm/aUZbZhl+op034/27RLy6ZtN20CRmCP6YAZKD2fjBwo9qr1SDj/fnj0m/Dgl+Ebx8GrL4FXXACWnNC6O9fVDCay7BjNsHMkzc7RDH1jXukdy9A/lqV/PEsiW5j28fURm8ZYkIZIgM6mGMesCtAQtakLB6iL2NRFvLo27JWakE3INuRzXEUk4c9SczxIfcTmub6pe+2P9DzC5b+7nF3pXbzt4LexpnYN9cF66kJ1xO34lA+Dq11c7eJoB1e7mLf9nPDGO0j8wzl8/MxXU9AFHO3guE6pLuhCab7gFiZq1yHv5qf0L84XlxfcAgVd8Gq/FJdnChnybr60PO/k9+ibd/Olx8w3hSptABQ3AgJmgIARKLUVS9AITkybQYKmNx8yQ6W2gBkgZIWmtIfMEEErSNgME7S8x4WtMEEziGXIW78qmBa88kI4/K3wP5d6N+J58ofe+Pydry53dAvKcTU7R9O8MJTmheEU3cNpeobT9Iyk6BlJ0zuaIe9M3RsPmAbLaoIsrwnxktYaTjo4yLKaIMviIZrjQZpiAZpjQRqiASxTjhgudvKtN0tKKQ5uifNc7xgAqXyKrz7xVX78/I/prO3khtfcwOFNh8/6+ZKPPcb27/6U+BvfwKEXfL6it5Jd7ZY2AIobBnnH2xjIubmJZW6enOPN55ycN+/myDv5Kcunq7NOdkpbcX40P0rOzZXaiu05J0fOze33OlmGVdoQCJkhwnaYsBkmZIUIW5Nqf1nEihC2woStMBE7QsSKTNT+dNSOErEi2KY9j/99AUBdB5z7A3j+1/DLi73b7r7kzXDqVdCwptzRzZt0zmHbUJKtgym27UqyfShVKj3D6SmH15WClniI9vowR3XU0/ayMG11IVprw7TWhmitDdEQDVT0d4tYWJLw5+CQ5XF++scetNbc9ORN3PH8Hbz3sPfy8aM+TsgKzfp5ptwU56rKvymOoYzSHjYVlMtc7U7ZEMg6WbIFr844GbIFv3ayZAqZUp9MIUPGyZApZEgX0qXp4nx/qp90IT2l5N38rOOyDZuoHfU2AOwIMTtWqmN2jKgdJRaIlebjgTixgFfH7Xhp3jYq6J9dKQ4+DdacCL//Bjx0rXdS3/EfhRM/DaGackc3K3nHZftQii0DSbYMJtk8mGTLYIKtgyl6xzJT+tZFbFY1RDhiRR1vOqKVjvoIK+ojdDSEaa0NE7Bkr1zMniT8OThkeZxEtsCO0QxPDTzFkc1HcvGxF8/pOXShwI5//DRuIsHKW76LGYseoGgXP0MZhKzQnDa29lfBLZAupEnlU15dSJHKp7y6OJ1PkcwnSRaSE9N+GcmM0DPeQzKfJJFPkC6k9/k3w1aYeCBOTaBmogS9ujZY65VA7cR0sJa6YB0xO1bxG5Evih32Bux5+bu9Q/yPXA9/vg1OvhSOPg8q5AjLaCpP10CCTQMJNvUn2DSQZPNAgu1DqSl76vURmzVNUU44qJE1jVFWNUVZ0xhlZWOE2nBlrItYHCThz0HxxL2/7Ryla6SL01afNufnGLj+BlKPP07b1V+Rm+JUEcuwvD3wQHxenq/gFkrJP5FLMJ4bZzw3TiKfYCw3RiLn1eO5ccZyY4zlxtiZ3Mlzw88xlhsjmU/OHKuyqA3WUh+qpy5YR32onvpgvVeH6mkINZRKY7iR2kAtpmHOy3otqJpW72Y8x30QfvU5+MU/wu9vhNf+P+83/wXY6NFa0z+epas/wca+cboGEnT1J+jqTzKYyJb6BUyDVY0RDm6J84aXLqezOUZnc5TOpih1ETkBUSwMSfhzsM5P+H/qeYGx3Bjr6tbN6fHj9z/Aru98h7pzzqH2zDMPRIiiSliGVdor3x95N894bpzR7Cij2VFGsiOMZEemTI9kRhjODrNpZBPDmWFGsiN4Y3xNZSiD+mA9jeFGGkONNIWbSqU50uzV4WaWRZYRsSMvdtXnX/sx8P57YOOv4b4r4M73w+9ugNd8Dg46ZV4Sv9aaHaMZL6n3J9jYl2Bj/zgb+xOMZybOeo+HLA5aFuM1hzRz0LIYBy2LsbY5xor6sJwUJ8pOEv4c1IZtWmtDPNX/VwDW1c8+4ee6u9lx6aWEDjuMls/ITXHEi2MbdmkvfbYc12E0N8pwZpihzBC7MrvYlfbKUGbIm87sYuvYVgbTg9OetxC1o6Xk3xJpYVlkmTcdbWF5dDktkRYaQg0LP0aEUnDw6+Gg13m33H3gS/CDs2HFsd6h/rWzS/yuq+kZSXvJvC/B830Juvq9JJ/MTVyt0hgNsK4lxpkvb2Pdsjjr/OTeHA8u7p9TRFWThD9HhyyP0zW2CcJwUN1Bs3qMm83Sc+FFALRffx1GcH5H6hNiNkzDLG0krGXtXvtqrRnLjTGYHmQgPcBAaqBU96X6GEgNsKFvA/2pfgp66nXdtmHTEmmhNdZKa9QrbbE2WqOttMfaWR5d7p0AeiAYJrz8nfDSs+HJ2+C3/+rdlGfFcXDSp2HdaaAUBcdl21CKjX3eb+wb+7y99U0DCTL5iZHhlsWDrGuJ8fb1HaxribFuWZyDlsVoiMpheFF9JOHP0SEtcR7r2k5LfRP1odmN8d335S+TefZZVnzzRgIdB2YsfiHmk1Kq9JPD2rqZNw5c7TKUGaIv1UdvspfeZC99yT56U73sTOzkDzv/wEB6AFdPJFGFojncTHu8nbZYG+2xdlbEVrAivoKOeAfN4eYXf06BFYBj3kfqsHcw/Mi/U7/h60RuewfdgU7+03wLt44dTcaZOArRVhvioJY4x3c2lvbW1y2LUxuRk+bE4iEJf44ObonD9l7ao7O79nf07rsZ+dHtNH7wfOKvfe0Bjk6IhWUoo/R7/+GN049DkXfz9Kf62ZHY4ZXkDnrGe+hJ9PCnvj/xyy2/nLJBYBs27bF2OuIdrKxZSUe8g454B6tqVtEWa9vjckWtNb1jGTb7Z8FvGkiyaSDB5oEkPSNpYDUWV3OG+Xs+bvyCy3Jf46OxVratfQ/GMe9hzYq2BR3PXIhykXf5HK1riWAE+4gbx+2zb3bjRnZ+/goi69fTfNFFCxCdEJWnmMDbY9Pf8jnv5ulN9PJC4gW6x7vpTnTTPd7NC+MvsKFvA6lCqtTXwKQu0EJYtUC+mWSinl3DtaRTDehCDWAQDZh0Nsc4dnU95zZ3cNCyGJ3NMVY3vYmg8SXY+CvqHr6Oumevga5vwsvfBX/3YWjc+88cQlQ7SfhzFI6Moow8Or98r/2cRJLuCy/CiEZp+9pXUZb8q4WYjm3YdNR00B5bQX80y/ahFNt0krZsihaVZEuil+1j20nThxHYRTYwiBnowwg8A5EcZgRiQMAI0RFfybr6TlbXrmZNzRrW1K5hdW0dYSs88QcPeaNXdvwJHv0WPHELPHYzrH0NHPN+b1mFXMsvxHySLDRH2xObARgbbZyxj9aa3sv/H7lt21h5yy3Yy5YtVHhCVCytNSOpPN3DabqHU7wwnCqNC//CUIoXhtNTbqVqGooV9WFWNjRydPtKVjdFWdMUYXVjlBX1EWxT0Z/qZ/v4draMbmHr2Fa2jm7l6cGn+fW2X0/5maAt2saaWm8DoLOuk87aTtY2rKXurG/DqVfChn+HP34PfvweiC6Do97tFdnrF4uIJPw52ji8EYAd/TNfPz38/R8wds8vaf7kJ4ke/3cLFZoQZZUruPSNZdjh36mtZyRNz0iaHSPFm7ykSeWm3oipJmTR0RBh3bI4r3tJCx0NEVb6pb0+jL2Pa9dboi20RFs4dvmxU9qzTpbtY96GwJbRLWwe3cyW0S1s6NtAxpkYvrYh1OAl/7q1dL7hs6xNDrN244M0PnId6uGveWf3H3kOHH4WRGZ/CaQQlUgS/hx1jXQRN1vYtqtAJu8QsqeeTZz+85/pu+YaYiefTOMHzy9TlELMH601Y5lC6VarfWNZ75aro95tV70kn5kyslxRfcSmrS7M6qYor1rXxIr6CO11YVbUh+loOHBDxwbNIOvq1+0xVoarXXYmd7J5ZDObR/0yspl7ttzDeG681K9m3UtYa0boHB9k7UNXsPbBz7O2/RUsO+xtqENPr5px+4WYTGm958hbi8X69ev1E088Ma/P+Za73kKQZfzh0TP5+QWv4mUrJvb0C0NDbDnrbJRlseYnd2LW7t8oakIshHTOYTCR9UuOgXFvemDcK/3jGQYSWfrHsmQnHWovqovYLK8J0VIT8u/O5t+lrc6bb6sLEwlUxz6F1prB9CCbRjexeWQzm0Y2sXnUq4ezw6V+MdelM1+gM7yMtS1H09l5KmuXH0VrtHXhBxsSYhpKqQ1a6/XTLauOT2OFyDk5to5t5azOE/kD8FzfeCnha8dhx8WX4AwNseqHt0myFwvKcTWj6TzDqRwjqRzDyTxDqRzDyRxDyRy7ptRZdiVyexxeL6oN2yyLB2mOBzlmZT3N8SAtfmJvqQmxzJ8PB6pw/P0ZKKVojjTTHGnm+Nbjpywbygx5GwDDm9jU83s2DzzJQ5lB7tp5P+y8H4Cwslgd76Cz8bDSOQJratewMr5SbpcsKoYk/DnYMroFRzsc3Xoo37cMnu+bOAQ4eOONJB95hOVXfoHw4dNfjyzE3hQcl0S2wHimwGg6z1gmz1i6wFg6z+ikMlKcTuUYSecZSXl9ZzpYF7AMGqMBGvyypjFCYyxIUyxIYyxAUyxAcyxEUzxAYzQot1zdTUOogYblDd55Ai8512vUmtHtD7P5rz9jU/cjbEr1sjk5zobhzfzCnBha11QmHfGO0gmDq2tWl+q6UF2Z1kgsVZLw56BrpAuAQxoO5qDmXv7W6yX8sXvvZfCbN1F71lnUvf3t5QxRLDCtNem8QzLrkMoVSGQLJLMOyWxx2qsT2QKJTKGU0MezBcYzeRIZfz6TnzJW+3QMBTVhm9qwTV0kQF0kwOqmaGm+PmJTHwlQF7FpiAaojwRojAUI26aM7z7flKJ21YkctepEjgIY2wGb7odN95Pa8iBb8uNsti221DSz1UqwOfsXHu55iLw7MQxxXbCOVTWrWFWzitU1q1lZs5KV8ZWsrFlJ1JbbZov5Jwl/DjYOb8QyLFbXrObQ1iT3PtPHo/c/Tv0/XUroiCNY/vnL5Yu1wjiuJpN3vFJwSeec0nw675DOeXWmNO2SzhVI5RxSflvKn0/nHJI5h3SuQDLnkMoWSOWdGfesdxcJmMSCFrGgRTxkEQ/ZtMRDxEMWNWHbq0NeXRu2qQnbE/MRm1jAwjDk/VWRatpKl/JFXJfDe5/i8K0Pw7ZHYNvvIDNCAdgRb2Zry8FsqWlhWyDANjfNozsf5e5Nd095usZQY2mEwY54R2nY4fZYO03hJvmeEftFEv4cdI10sbpmNbZp87GT1/Lcxh0kLv4UASNA21evlZviTOK6mpzjknNc8gWvzhVc8o5LtlCc1uQKLjnHIZv3+hSXFWtv2im1ZQuOV+cnpjN5v82fzuRdsnmHTMEh7+zfSamRgEnYNgkHTCIBk0jAIhIwqYvYhAMWsaDXFg2Ypflo0CIS8BJ6NGj6tVdiQQtTkvXSYBjQ9nKvnHABuC70P4vV/RgruzewsvtxTup6ZKJ/vI3U8peyvXEl2yO1bLctXnBSbE/08FjvY/x888+nPH3QDNIWa/NKtK10Y6LWaCvLo8tpjjTvMfywECAJf066Rro4oukIANY2Rrhxy10kMiNccsKHSd/+PJe/2WJtU4z6qE0saM16K1xrjauh4Lo4rqbgahxH42g9Zb64PO8U210Krqbgz+dd15/2kmnBn/f6TPQt9is4Lnl/2eT+eUeTd1wKrkuuMNGec9xS33wxmTsu+cJu834888UyFEHLIOCXkG0S9OuAaRAJWDREDYKWSdCeujxkmYT8tmIdtk2vDnjLwwG/lJYZsgcl5o9hwPKXemX933tt6RHofQp2PgU7nyTS+xSHdv2GQ7X/s44yoaETmtaRXX4sPbEmuoNhegxNTyFJT3InPYkenhl8hpHsyJQ/V7w5UUt04vbFxdsZF+970BxupjZYK+/zJUYS/iwl80l6Ej2cve5s3GyW3i9cSfrhh2i/4go+d+wpXHzHU7z/3x8v9Q+YBvGQhVKUDvlqvBOzXO0dai4m9PlMjnNlmwrLMLBNhW0aWJPmLdPAMhQBy/CWGYpIwPLnvf62aRAwDWzLe1zA8uf9toBplB7v9fPqYLHNmnhM0J5YFrTM0jLZMxaLTrgO1pzklaJCFgY3Qv9fYeCvMPAc7NpEsOs3dDpZOov9lAm1K6BuJdQeTqqxhd5wjJ12gF7l0qtz9ObH6Uv1s3V0K4/tfIzx/PgeIViGRUOogcZQIw1hvw41UBesoyHUQH2onrpgXemuiTWBGixDUkY1q6pXTyn1BuB6wAT+TWv9lYX626UT9jL1bHvXu8k8/TRNH/so9eeewwnArz55Eo9t2cWuRI7hVI6hZJ7xTH5S7F5tGQaGUpgGGIbCMhSmYfi1KtVmad6Y0uYlZK+9OG8aXvI1DYVtGBN9/CRdTOLF6WLyNg0lW/hCVAorOHEkYDLXgbEe2LUJRrbB8DYY2e6VLQ8RGd9Jp3YmNggAUBBphGgzRJtIhTsZCMcYsIIMWiaDSjOAw5CbZVchxVCyn03DGxnOjpJ19hxAqShmx6gJ1FATrKEmUEM8ECdmx7w6ECNmx4jaUWJ2jIgdIWpHiVgRInaEiBUhbIUJW+EXf/tjsV+qJuErpUzgRuBUoBt4XCl1t9b62YX4+13DXRy5yaX5m9eQc2HFjd8gfsoppeWxoMVrD21ZiFCEEEuJYXp783Urp1/uFCDZ710pMN4Lib6JkhyE1C4i/X9jVWqQVekRvGON09NAOhBhOBRnOBhhxA4yagUYsSxGDYMxDMbyDuP5IcYYYJsukHALJHSehJub9SrZhkXYDBE2Q4TMIEErSMifDvjTATNI0PJq2wwQKBYjgG3Y2Kbt1YaNZVjevPKmTcPEMiwsZXnzyiy1m8ov/rShDCzDwlBGab5YTOVd4WIqE0X17yBVTcIHjgO6tNabAZRSPwLOBBYk4atb7+Syn7gED2llxQ03EFi1aiH+rBBC7J1peVcJ1LTtu6/rQnYU0sNeyYxBZhSyXq2yCSK5BJHsOO25BORSkE9CNgX5FOTTXin4tTOR5B0grRQJwyBl+LVSpAyDpFKkDYO0UqQMRVoZpA1FRimyyqszhjed9NuKJa8UOb/ky5xwlQYDUHh1cVpNmVZT2vcsasr8CZEOPv+O/16Q+Ksp4bcDL0ya7wb2uDONUupDwIcAVq6cYYt4P4RqGthywkpOv/FHGOHwvh8ghBCVxjAgXO+V+eC6XtIvZDALGWJOjlgh57U5We/og1Ocz4NbANevnYJXa8dvd7w+2vGmteM9v3a9ae2inQIF7ZBz8+S1Q8715vOu49W6QEG7fnEoaBcXl4LWXhsurtY4uN45VH7toilojevPa8DBay8Wrb02DWg0LuD6fV2m1pP7aP+x7LHcc1B44e6mWk0Jf1a01jcDN4M3lv58Pe8bL/kGQNUf0hFCiHljGGCEwA4tyJ9TgO0XMXfVlPB7gI5J8yv8tgUhiV4IIUQ1q6ZBsx8H1iml1iilAsC5wN37eIwQQgghqKI9fK11QSl1AfArvMvybtFaP1PmsIQQQoiqUDUJH0BrfQ9wT7njEEIIIapNNR3SF0IIIcR+koQvhBBCLAGS8IUQQoglQBK+EEIIsQRIwhdCCCGWAEn4QgghxBIgCV8IIYRYAiThCyGEEEuAJHwhhBBiCVBaz9sN5SqOUmoA2DaPT9kEDM7j85XTYlmXxbIeIOtSqRbLuiyW9QBZl71ZpbVunm7Bok74800p9YTWen2545gPi2VdFst6gKxLpVos67JY1gNkXfaXHNIXQgghlgBJ+EIIIcQSIAl/bm4udwDzaLGsy2JZD5B1qVSLZV0Wy3qArMt+kd/whRBCiCVA9vCFEEKIJUAS/jSUUm9QSj2nlOpSSl06zfKgUup2f/kflFKrFz7KvVNKdSilHlBKPauUekYpddE0fU5WSo0qpf7sl8vLEetsKKW2KqX+4sf5xDTLlVLqBv81eUopdXQ54twXpdQhk/7ff1ZKjSmlPrFbn4p9XZRStyil+pVST09qa1BK3auU2ujX9TM89jy/z0al1HkLF/X0ZliXf1FK/c1/D/1MKVU3w2P3+n5cSDOsxxVKqZ5J76HTZ3jsXr/rFtoM63L7pPXYqpT68wyPrZjXxI9n2u/gsn5etNZSJhXABDYBnUAAeBI4bLc+HwO+5U+fC9xe7rinWY9W4Gh/Og48P816nAz8d7ljneX6bAWa9rL8dOCXgAKOB/5Q7phnsU4m0It33WxVvC7AScDRwNOT2q4BLvWnLwWunuZxDcBmv673p+srcF1OAyx/+urp1sVfttf3YwWsxxXAp/fxuH1+11XCuuy2/KvA5ZX+mvjxTPsdXM7Pi+zh7+k4oEtrvVlrnQN+BJy5W58zgVv96TuBU5RSagFj3Cet9U6t9R/96XHgr0B7eaM6oM7CUip/AAAGsUlEQVQEvqc9jwJ1SqnWcge1D6cAm7TW8zk41AGltf4tMLRb8+TPw63AW6Z56OuBe7XWQ1rrYeBe4A0HLNBZmG5dtNa/1loX/NlHgRULHtgczfCazMZsvusW1N7Wxf+OfQfwwwUNaj/t5Tu4bJ8XSfh7agdemDTfzZ6JstTH/3IYBRoXJLr94P/kcBTwh2kWv0Ip9aRS6pdKqcMXNLC50cCvlVIblFIfmmb5bF63SnMuM395VcvrAtCitd7pT/cCLdP0qcbX5+/xjhpNZ1/vx0pwgf/TxC0zHDauttfkRKBPa71xhuUV+5rs9h1cts+LJPxFTikVA34CfEJrPbbb4j/iHU4+Evg6cNdCxzcHr9JaHw28EfgHpdRJ5Q7oxVBKBYAzgDumWVxNr8sU2jseWfWX/iilPgsUgB/M0KXS3483AWuBlwM78Q6FV7t3sve9+4p8Tfb2HbzQnxdJ+HvqATomza/w26bto5SygFpg14JENwdKKRvvjfYDrfVPd1+utR7TWif86XsAWynVtMBhzorWusev+4Gf4R2OnGw2r1sleSPwR6113+4Lqul18fUVfz7x6/5p+lTN66OUeh/wJuBd/hfyHmbxfiwrrXWf1trRWrvAd5g+vmp6TSzgLOD2mfpU4msyw3dw2T4vkvD39DiwTim1xt8LOxe4e7c+dwPFsybfBtw/0xdDufi/d30X+KvW+msz9FlePPdAKXUc3vuhEjdcokqpeHEa78Sqp3frdjfwXuU5HhiddNisEs24t1Itr8skkz8P5wH/NU2fXwGnKaXq/cPLp/ltFUUp9QbgEuAMrXVqhj6zeT+W1W7nr7yV6eObzXddpXgd8Detdfd0CyvxNdnLd3D5Pi/lPpOxEgveGd/P453B+lm/7Uq8LwGAEN6h2C7gMaCz3DFPsw6vwjtU9BTwZ7+cDnwE+Ijf5wLgGbyzcx8FTih33DOsS6cf45N+vMXXZPK6KOBG/zX7C7C+3HHvZX2ieAm8dlJbVbwueBspO4E83u+KH8A7f+U3wEbgPqDB77se+LdJj/17/zPTBby/QtelC++30+Jnpng1Thtwz97ejxW2Hv/pfw6ewkswrbuvhz+/x3ddpa2L3/4fxc/HpL4V+5r4Mc30HVy2z4uMtCeEEEIsAXJIXwghhFgCJOELIYQQS4AkfCGEEGIJkIQvhBBCLAGS8IUQQoglQBK+EGKvlFIPKqW0PyCNEKJKScIXYpHzk/X+lAfLHbsQYv5Y5Q5ACHHA7TF8r68BsIEM3g2gdle8a9l24LkZ+gghqoQMvCPEEuXvwb8auFVr/b7yRiOEONDkkL4QQgixBEjCF0Ls1Uwn7SmlTvbbt/rzr1dK3aeUGlJKjSil7lVKvWJS/1ql1D8rpZ5XSqWVUi8opa5WSoX38fdfpZT6kVKqWymVVUrt8v/OO4s3GRJC7Jv8hi+EeNGUUh8DvoF3s5AEUIN3h7NXKaVeh3eDlvuBlwJJvJ2NFXh3pjsc73a00z3v1X6fojGgHjjFL2copd6lvdvACiH2QvbwhRAvVjNwLfBloFFrXQusAX6Pd2fJa4Fv4Z0geCIQ98v5QAH4P0qp03d/UqXURXjJvg/4EFDnP3cU71auvX79Twdy5YRYLCThCyFerAhwm9b6s1rrEQCt9Va8ZKyBY4E3A2/SWj+sPTmt9XeB7/nP8bbJT6iUqgO+iHcFweu11t/RWo/6z53WWt8OnOU//8X+/dyFEHshCV8IMR++vHuD1no73j2/Ae7QWndN87jf+PVLd2s/G4gB92mtn5zuD2qtfw9swTvEf8z+BC3EUiK/4QshXqwME4l9d/3AwcDTMywvjhFQv1v7CX79WqVU717+doNfd+D9hCCEmIEkfCHEi9WnZx7Qw/HrnftYbu/W3urXEb/sy2z6CLGkScIXQlSi4s+N12utP1HWSIRYJOQ3fCFEJSoe6l9Z1iiEWEQk4QshKlHx9/iT9zUwjxBidiThCyEq0R14A/TUA5fvraNSavcT/oQQ05CEL4SoOFrrXcBl/uylSqnvKKUOLi5XSoWVUicqpW4CfleWIIWoMnLSnhCiImmtv66UqgWuxBuV73ylVBLIAbVM7LBsLU+EQlQX2cMXQlQsrfUXgSOBm/Gu9TfwhtbdCfwKb+jdE8sWoBBVRM18+awQQgghFgvZwxdCCCGWAEn4QgghxBIgCV8IIYRYAiThCyGEEEuAJHwhhBBiCZCEL4QQQiwBkvCFEEKIJUASvhBCCLEESMIXQgghlgBJ+EIIIcQSIAlfCCGEWAL+P+oOkWZptOhdAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import pylab as plt\n",
    "from scipy.integrate import odeint\n",
    "\n",
    "# initial conditions\n",
    "y0 = (10, 1, 0, 0) # (S0, E0, ES0, P0)\n",
    "kf, kb, kc = (1., 1e-2, 1.)\n",
    "\n",
    "time = np.linspace(1e-2, 20, 101)\n",
    "res = odeint(conversion, y0, time, args=(kf, kb, kc))\n",
    "dS, dE, dES, dP = zip(*res)\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))\n",
    "ax.plot(time, dE, label='$dE$')\n",
    "ax.plot(time, dS, label='$dS$')\n",
    "ax.plot(time, dES, label='$dES$')\n",
    "ax.plot(time, dP, label='$dP$')\n",
    "ax.set_xlabel('Time', fontsize=24)\n",
    "ax.set_ylabel('Concentration', fontsize=24)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "c2_bqXG6L2Hu"
   },
   "source": [
    "#### Extra - C++ version"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ybsPCd89L7CT"
   },
   "source": [
    "We propose the C++ version of the same model.\n",
    "\n",
    "In the version proposed above we have used the *odeint* function provided by *scipy*, while in C++ we need to implement the integration by hand.\n",
    "In all the previous examples we have used the Euler method for the integration of ODEs; in this case we introduce a more advanced technique of ODE integration: the **Runge-Kutta** (RK) method.\n",
    "\n",
    "The Runge-Kutta method is a mathematical algorithm used to solve systems of ordinary differential equations (ODEs). \n",
    "The general form of these equations is as follows:\n",
    "\n",
    "$$\n",
    "\\dot{x} = f(t, x) \\quad \\quad x(t_0) = x_0\n",
    "$$\n",
    "\n",
    "Where $x$ is either a scalar or vector. \n",
    "The 4th order Runge-Kutta method (yes, there is different level of precision of the same method and the RK4 is the most common one in numerical simulation due to its good balance between numerical efficiency and numerical stability) is given by:\n",
    "\n",
    "$$\n",
    "x_{i+1} = x_i + (k_1 + 2(k_2 + k_3) + k_4)/6 \\quad \\quad t_{i+1} = t_i + h\n",
    "$$\n",
    "\n",
    "where $h > 0$ is a step size parameter, $i=1, 2, 3, \\cdots$ and:\n",
    "\n",
    "$$\n",
    "\\begin{cases}\n",
    "    k_1 = f(t_i, x_i)h                               \\\\\n",
    "    k_2 = f(t_i + \\frac{h}{2}, x_i + \\frac{k_1}{2})h \\\\\n",
    "    k_3 = f(t_i + \\frac{h}{2}, x_i + \\frac{k_2}{2})h \\\\\n",
    "    k_4 = f(t_i + h, x_i + k_3)h                     \\\\\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "The Runge-Kutta method offers greater accuracy than the method of multiplying each function in the ODEs by a step size parameter and adding the results to the current values in $x$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3JGIINTHMAA4"
   },
   "source": [
    "In the MM model we have a system of 4 ODEs that we have to integrate simultaneously.\n",
    "Therefore we need to re-apply the same formalism on each equation.\n",
    "\n",
    "Lets see how our C++ code (could) appear (ref. previous Euler implementation)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "_-wRqmu-MCr1"
   },
   "source": [
    "```c++\n",
    "// g++ michelis_menten_rk4.cpp -std=c++14 -O3 -o michelis_menten_rk4\n",
    "\n",
    "#include <memory>\n",
    "#include <algorithm>\n",
    "#include <cmath>\n",
    "#include <type_traits>\n",
    "#include <cassert>\n",
    "\n",
    "template < class type >\n",
    "using is_floating_point = typename std :: enable_if < std :: is_floating_point < type > :: value > :: type *;\n",
    "\n",
    "template < class type, is_floating_point < type > = nullptr >\n",
    "using array = std :: unique_ptr < type[] >;\n",
    "\n",
    "\n",
    "/**\n",
    "* @brief Michaelis Menten kinetic\n",
    "*\n",
    "* @param x List of time points.\n",
    "* @param s0 Initial condition of the substrate.\n",
    "* @param e0 Initial condition of the enzyme.\n",
    "* @param es0 Initial condition of the substrate+enzyme.\n",
    "* @param p0 Initial condition of the product.\n",
    "* @param kf Constant of the forward reaction.\n",
    "* @param kb Constant of the backward reaction.\n",
    "* @param kc Constant of the product reaction.\n",
    "* @param S The resulting substrate array.\n",
    "* @param E The resulting enzyme array.\n",
    "* @param ES The resulting enzyme+substrate array.\n",
    "* @param P The resulting product array.\n",
    "*\n",
    "* @tparam type Data-type of arrays\n",
    "* @tparam Length of time points.\n",
    "*\n",
    "*/\n",
    "template < class type, int32_t N >\n",
    "void MichaelisMenten (const array < type > & x,\n",
    "                      const type & s0, const type & e0, const type & es0, const type & p0,\n",
    "                      const type & kf, const type & kb, const type &kc,\n",
    "                      array < type > & S, array < type > & E, array < type > & ES, array < type > & P\n",
    "                     )\n",
    "{\n",
    "  // determine the interval as diff\n",
    "  const type dx = x[1] - x[0]; // Note: we are assuming it is constant!!\n",
    "\n",
    "  // Set the initial condition\n",
    "  S[0] = s0;\n",
    "  E[0] = e0;\n",
    "  ES[0] = es0;\n",
    "  P[0] = p0;\n",
    "\n",
    "  // set the equation functions\n",
    "  auto dS = [](const type & S, const type & E, const type & ES,\n",
    "               const type & kf, const type & kb)\n",
    "            {\n",
    "              return -kf * S * E + kb * ES;\n",
    "            };\n",
    "  auto dE = [](const type & S, const type & E, const type & ES,\n",
    "               const type & kf, const type & kb, const type & kc)\n",
    "            {\n",
    "              return -kf * S * E + kb * ES + kc * ES;\n",
    "            };\n",
    "  auto dES = [](const type & S, const type & E, const type & ES,\n",
    "                const type & kf, const type & kb, const type & kc)\n",
    "             {\n",
    "               return kf * S * E - kb * ES - kc * ES;\n",
    "             };\n",
    "  auto dP = [](const type & ES, const type & kc)\n",
    "            {\n",
    "              return kc * ES;\n",
    "            };\n",
    "\n",
    "  // Integrate the equation using the RK4 method\n",
    "  for (int32_t i = 0; i < N - 1; ++i)\n",
    "  {\n",
    "    const type ks1  = dx * dS( S[i], E[i], ES[i], kf, kb);\n",
    "    const type ke1  = dx * dE( S[i], E[i], ES[i], kf, kb, kc);\n",
    "    const type kes1 = dx * dES(S[i], E[i], ES[i], kf, kb, kc);\n",
    "    const type kp1  = dx * dP(ES[i], kc);\n",
    "\n",
    "    const type ks2  = dx * dS( S[i] + .5 * ks1, E[i] + .5 * ke1, ES[i] + .5 * kes1, kf, kb);\n",
    "    const type ke2  = dx * dE( S[i] + .5 * ks1, E[i] + .5 * ke1, ES[i] + .5 * kes1, kf, kb, kc);\n",
    "    const type kes2 = dx * dES(S[i] + .5 * ks1, E[i] + .5 * ke1, ES[i] + .5 * kes1, kf, kb, kc);\n",
    "    const type kp2  = dx * dP(ES[i] + .5 * kes1, kc);\n",
    "\n",
    "    const type ks3  = dx * dS( S[i] + .5 * ks2, E[i] + .5 * ke2, ES[i] + .5 * kes2, kf, kb);\n",
    "    const type ke3  = dx * dE( S[i] + .5 * ks2, E[i] + .5 * ke2, ES[i] + .5 * kes2, kf, kb, kc);\n",
    "    const type kes3 = dx * dES(S[i] + .5 * ks2, E[i] + .5 * ke2, ES[i] + .5 * kes2, kf, kb, kc);\n",
    "    const type kp3  = dx * dP(ES[i] + .5 * kes2, kc);\n",
    "\n",
    "    const type ks4  = dx * dS( S[i] + ks3, E[i] + ke3, ES[i] + kes3, kf, kb);\n",
    "    const type ke4  = dx * dE( S[i] + ks3, E[i] + ke3, ES[i] + kes3, kf, kb, kc);\n",
    "    const type kes4 = dx * dES(S[i] + ks3, E[i] + ke3, ES[i] + kes3, kf, kb, kc);\n",
    "    const type kp4  = dx * dP(ES[i] + kes3, kc);\n",
    "\n",
    "    S[i + 1]  = S[i]  + type(1. / 6.) * (ks1  + type(2.) * ks2  + type(2.) * ks3  + ks4);\n",
    "    E[i + 1]  = E[i]  + type(1. / 6.) * (ke1  + type(2.) * ke2  + type(2.) * ke3  + ke4);\n",
    "    ES[i + 1] = ES[i] + type(1. / 6.) * (kes1 + type(2.) * kes2 + type(2.) * kes3 + kes4);\n",
    "    P[i + 1]  = P[i]  + type(1. / 6.) * (kp1  + type(2.) * kp2  + type(2.) * kp3  + kp4);\n",
    "  }\n",
    "}\n",
    "\n",
    "\n",
    "int32_t main (/*int32_t argc, char ** argv*/)\n",
    "{\n",
    "  const float s0  = 10.f;\n",
    "  const float e0  = 1.f;\n",
    "  const float es0 = 0.f;\n",
    "  const float p0  = 0.f;\n",
    "\n",
    "  const float kf = 1.f;\n",
    "  const float kb = 1e-2f;\n",
    "  const float kc = 1.f;\n",
    "\n",
    "  const int32_t iterations = 101;\n",
    "  const float dt = (20.f - 1e-2f) / iterations;\n",
    "\n",
    "  array < float > time(new float[iterations]);\n",
    "  std :: generate_n(time.get(), iterations, [n = 0, dt] () mutable { return 1e-2f + dt * n++; });\n",
    "\n",
    "  array < float > S(new float[iterations]);\n",
    "  array < float > E(new float[iterations]);\n",
    "  array < float > ES(new float[iterations]);\n",
    "  array < float > P(new float[iterations]);\n",
    "\n",
    "  MichaelisMenten < float, iterations >(time, s0, e0, es0, p0,\n",
    "                                        kf, kb, kc,\n",
    "                                        S, E, ES, P);\n",
    "\n",
    "  return 0;\n",
    "}\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "e6QVRNvDMFYD"
   },
   "source": [
    "The code is certainly longer than the previous ones, since the RK4 method is more verbose and it requires multiple steps of integration.\n",
    "However the implementation is very simple also in C++ and its transposition in Python is left as exercise.\n",
    "\n",
    ">⚠️**Note:** All the numerical ODE integration require an iterative process by definition.\n",
    "> The iteration is always traslated into a *for* or *while* loop and its parallelization is very (**very!**) difficult, if not impossible!\n",
    "> Everithing involve a *for/while* loop is extremely slow in Python and also *scipy* use particular hacking under the hood to improve the numerical efficiency.\n",
    "> If you want to try to re-implement the RK4 method in Python, please pay attention to a good usage of NumPy library!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "eO4QEQR6MGnu"
   },
   "source": [
    "# Hill kinetic"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "UHZ_FpQKMI8M"
   },
   "source": [
    "In biochemistry and pharmacology, the Hill equation refers to two closely related equations that reflect the binding of ligands to macromolecules, as a function of the ligand concentration. \n",
    "A ligand is \"a substance that forms a complex with a biomolecule to serve a biological purpose\" (ligand definition), and a macromolecule is a very large molecule, such as a protein, with a complex structure of components (macromolecule definition). \n",
    "Protein-ligand binding typically changes the structure of the target protein, thereby changing its function in a cell. \n",
    "\n",
    "The Hill kinetic borns aiming to explain the binding curve between the 4 haemoglobin (Hb) sub-units and the oxygen (O).\n",
    "Each bound between these two substances tends to promote the bind of the other sub-units, with a change of Hb functionality.\n",
    "\n",
    "In general, the Michaelis-Menten kinetic produces *hyperbolic* trends, while there is a class of enzymes which describes *sigmoidal* curves.\n",
    "Sigmoidal profiles are due to enzymes with a *positive co-operative* behaviour.\n",
    "\n",
    "Cooperativity is a phenomenon displayed by systems involving identical or near-identical elements, which act dependently of each other, relative to a hypothetical standard non-interacting system in which the individual elements are acting independently.\n",
    "One manifestation of this is enzymes or receptors that have multiple binding sites where the **affinity** of the binding sites for a ligand is apparently increased, positive cooperativity, or decreased, negative cooperativity, upon the binding of a ligand to a binding site.\n",
    "A cooperative binding requires that the macro-molecula (Hb in this case) has more than one binding site (sub-units).\n",
    "\n",
    "The positive cooperation implies an **allosteric regutation**, i.e the regulation of an enzyme or a protein through a molecula called **effector**.\n",
    "In biochemistry, allosteric regulation (or allosteric control) is the regulation of an enzyme by binding an effector molecule at a site other than the enzyme's active site.\n",
    "The site to which the effector binds is termed the allosteric site or regulatory site. \n",
    "Allosteric sites allow effectors to bind to the protein, often resulting in a conformational change involving protein dynamics.\n",
    "Effectors that enhance the protein's activity are referred to as allosteric activators, whereas those that decrease the protein's activity are called allosteric inhibitors. \n",
    "\n",
    "The **cooperativity degree** is determined by the *Hill equation* in kinetic reactions which follows the Michaelis-Menten equations, i.e for the description of enzymes with more than one sub-units or for multi-enzymes systems.\n",
    "The Hill equation considers the allosteric regulation of different sites using using the formula\n",
    "\n",
    "$$\n",
    "\\theta = \\frac{[L]^n}{k_d + [L]^n} = \\frac{[L]^n}{k_a^n + [L]^n}\n",
    "$$\n",
    "\n",
    "where $\\theta$ is the fraction of the receptor protein concentration that is bound by the ligand; $[L]$ is the total ligand concentration; $k_d$ is the apparent dissociation constant derived from the law of mass action; $k_a$ is the ligand concentration producing half occupation; $n$ is the Hill coefficient which describe the cooperativity ($n > 1 \\rightarrow$ positive cooperation; $n = 1 \\rightarrow$ Michaelis-Menten equation).\n",
    "\n",
    "When $n = 1$ we have an equivalence between all the sites.\n",
    "In general, $n$ should be equal for all the sites and only fitting  real data we can obtain the true value of the Hill coefficient (e.g. in the haemoglobin case we can obtain a value $n \\sim 2, 3$).\n",
    "\n",
    "The Hill equation could be obtain by the system of kinetic equations\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    [E] + [S] \\underset{k_b}{\\stackrel{k_f}{\\rightleftharpoons}} [ES] \\\\\n",
    "    [ES] + [S] \\underset{k_b}{\\stackrel{k_f}{\\rightleftharpoons}} [ES]_2 \\\\\n",
    "    [ES]_2 + [S] \\underset{k_b}{\\stackrel{k_f}{\\rightleftharpoons}} [ES]_3 \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where $[ES]$ will follow a Hill kinetic.\n",
    "\n",
    "The $n$ and $k_i$ parameters are hard to deduce in chemistry and it is preferable to obtain them via experimental fitting.\n",
    "To this purpose is preferable to transform the Hill equation in its logarithmic form, i.e\n",
    "\n",
    "$$\n",
    "log\\left(\\frac{\\theta}{1 - \\theta}\\right) = n log[L] - \\log(k_d)\n",
    "$$\n",
    "\n",
    "in which $n$ is related to the sigmoidal form of the equation, while $k_d$ is related to the velocity of site saturation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ynfy4hJfMNB-"
   },
   "source": [
    "Starting from the \"classical\" enzymatic kinetic we had 4 objects to combine: $[E]$, $[S]$, $[ES]$, $[P]$.\n",
    "Now, the Hill hypothesis can be re-formulated as the introduction of an *activator* in the same kinetic equation, i.e\n",
    "\n",
    "$$\n",
    "[E] + [S] \\rightarrow [ES] + [S] \\rightarrow [ES] + [P]\n",
    "$$\n",
    "\n",
    "In this way the enzyme could be activated only when a given amount of substrate is already present in my system.\n",
    "This structure introduces a level of sensitivity about the amount of substrate and model a wide range of real situations.\n",
    "More step we insert in the middle of reaction and more sensitivity we add to the system.\n",
    "\n",
    "This situation is very close to the Michaelis-Menten one, but we have to include an exponential term to $[S]$ and $k_m$.\n",
    "This exponent is related to number of sites required to verify the transformation in products.\n",
    "\n",
    "Lets try to use Python code to simulate the different behaviors of Hill equation, varying its parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 409
    },
    "executionInfo": {
     "elapsed": 832,
     "status": "ok",
     "timestamp": 1632141111544,
     "user": {
      "displayName": "nico nico",
      "photoUrl": "https://lh3.googleusercontent.com/a/default-user=s64",
      "userId": "04264570794863623616"
     },
     "user_tz": -120
    },
    "id": "lUj4WU8uL-Nn",
    "outputId": "1e92a397-ecd7-4194-83ab-fc5a539d9533"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7eff4c67f190>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAewAAAF2CAYAAACs16tFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdeXxU1f3/8deZrGTfCcmQhS0sAQIEEFnEgtaF4la1ot9q1dp+W22t1m/110Vb22pb29rWrra1LtVa3BdEQUUQlR3ZdwIkARIC2cl+fn/MBANMNjKTTMj7+XjkkczcM+d+sJR3zr3nnmOstYiIiIh/c/R0ASIiItI+BbaIiEgvoMAWERHpBRTYIiIivYACW0REpBdQYIuIiPQCCmwRHzHGZBhjrDEmtyOvz/AcM919JHirbhHxTwpskU4yxvzLGPOGh/dz3eGZ4X7rADAAWN/Fc1n3V70xZo8x5hFjTPiZ9tldjDGXG2M+NsaUGmMqjTHbjDF/7+m6RHorBbaIj1hrG621h6y1DV3sajGu4B8E/AD4BvBIV+vzJWPMLGA+8DpwDjAOuAcwPVmXSG+mwBbxEW9c8nardQf/AWvts8C/gctPaTPWGLPCGFNtjFltjBnfoo54Y8xzxph8Y8xxY8xmY8xXTql1hjHmE/dIuMwYs9IYk93i+LnGmA/c/RcYY/5sjIlqo+YvACustT+31m6z1u601r5urb2lRZ83uc/3BWPMDmNMjTHmfWPMoBZtBhtjXjXGHDLGVBlj1hpj5pxSe7Ax5ufGmH3GmFr3VYhvtTg+0hjzpjGmwhhT5P5vkdzB//YifkOBLdL7HAeCTnnvIeBeYDxQAvzbGNM8mg0F1gJzgFHA74C/ukfBGGMCgVeBD4GxwGTgUaDRfXw08A7wmvv4lUAO8M82ajwEDDfGjG3nzxIC3A98BZgCBAAvtag9AngLuMB97hfdx4e36ONJ4MvAXcAI4Bag1F37AGApsAmYBMx29/mqMUb//knvYq3Vl7701Ykv4F9AA1B5ylc1YIEMd7sM9+vcjrxu41xvtHg9CTgCPO9+PdPdx+dbtJnqfs/ZRr//Af7u/jnO3f68Vto+BfzjlPdy3J9JauUz4cCb7jYHgBeArwMRLdrc5D4+tcV76bh+UZjdRu2fAD9w/zzU3cdFrbT9CfDuKe/Fuj8zqaf/LulLX5350m+YImdmKa7Qavk1z0fnush96bgG+Nh97jtOabOhxc+F7u9JAMaYAGPM940xG4wxJcaYSlyj5DQAa+1RXL8YvO2+dHyXMSatRX8TgBvcNVS6P7/cfWywp4KttVXW2kuBIcCPcY14HwI2G2P6t2jaBKxs8bl97vpHumsPN8b80hizxRhzzH3u3Obacd0bbwLeb+W/3QRgxim1H2irdhF/FdjTBYj0UtXW2l0t3zDGxPjoXEuB24B6oNBaW++hTcv3mrfga/6F/LvA3cC3gY24rgb8HHegA1hrv2KMeRS4CJgL/MwYc7m19m13P38HfuvhvAVtFW6t3Q3sBv5ujPkZsAP4X+ABD/V68oi7pu8CO3FdxXgKCG7rvC04cI30v+vh2OEO9iHiFxTYIv7vtF8OOmka8Lq19mkA9/3hYbjv8zaz1n4KfAr8whjzFnAj8Dau+9+julgDQB6uwI1o8Z4D12X+j9y1pQEpwNYWtT9lrX3RfTwU18h4h/v4encf5wMLPZxzLXANsK+VX3REeg1dEhc5++0AZhljprknaz0GZDYfNMZkGmMeds8ETzfGnA+MAba4m/wCmGSM+YsxZpwxZogxZo4x5q+tndAY84D7UvZMd//jcE1Si8A1ea1ZA/CoMWaKMSYH1wSyzbgeZWuu/QpjzHj35LdncE2iA8BauwP4L64R/FXuc003xvyPu8kfgWjgeWPMZGPMIGPMbGPM34wxkWfyH1OkpyiwRc5+P8V1n/gtXJfXq3A9GtasGteIez6ugHzSffwXANbaDcAMXJPkPsA1Cn+Iti8pf4Drl4IncY2W33Z/fq61dmmLdrXAz3Bd5l6B69+kK621zZfJ7wKKgGXu+j9x/9zSl4Fngd8D23Ddj492116IaxJeE64R+GZcIV7r/hLpNcxn/78QEek+xpibgMestRHttRURjbBFRER6BQW2iIhIL6BL4iIiIr2ARtgiIiK9gAJbRESkF/DbhVMSEhJsRkZGT5chIiLSbdasWXPEWpvo6ZjfBnZGRgarV6/u6TJERES6jTFmX2vHdElcRESkF1Bgi4iI9AIKbBERkV7Ab+9hi4icTerr68nPz6empqanSxE/EBoaitPpJCgoqMOfUWCLiHSD/Px8IiMjycjIwLXDqfRV1lpKSkrIz88nMzOz/Q+46ZK4iEg3qKmpIT4+XmEtGGOIj4/v9NUWBbaISDdRWEuzM/m7oMAWERHpBRTYIiIivYBXAtsYc5ExZrsxZpcx5l4Px0OMMc+7j68wxmR447wiItIzbr75ZpKSksjOzj7jPowx3HDDDSdeNzQ0kJiYyJw5czp0vCsiIiLabVNaWsqf/vSnE6/PPffcLp+3K7oc2MaYAOCPwMXASOA6Y8zIU5rdAhyz1g4Bfgv8oqvnFRGRnnPTTTexcOHCLvURHh7Opk2bOH78OACLFi0iNTW1w8d97dTA/uijj7rt3J54Y4Q9Cdhlrd1jra0D/gNcdkqby4An3T+/AMwy3Tj7oq6hkX+/9TNe/eRfrCj8hNWHVrO+aD0bizeyt2wv1fXV3VWKiEiP2rx5M7Nnz2bYsGE8+OCD3HHHHaxatarT/cyYMYO4uLgu13PJJZfw5ptvAvDcc89x3XXXdep4s6qqKi699FLGjh1LdnY2zz//PAC/+c1vyM7OJjs7m0cfffS0z+Xl5Z10leCRRx7hgQceAODee+9l9+7d5OTkcM8995w0KvfUb15eHiNGjOCrX/0qo0aN4sILLzzxy4Y3eOM57FTgQIvX+cDk1tpYaxuMMWVAPHCkZSNjzG3AbQBpaWleKM2lobaSh4v+A0XAds9tksKSGJ80ntnpszl/4PkEBwR77fwiIi39+PXNbCks92qfI1OiuP8Lo9psU1NTw9VXX838+fMZNGgQw4cPZ8KECUycOPFEm+nTp1NRUXHaZx955BFmz57t1ZoBvvSlL/GTn/yEOXPmsGHDBm6++WaWLVvW4ePNFi5cSEpKyolwLysrY82aNTzxxBOsWLECay2TJ0/mvPPOY9y4cR2q7eGHH2bTpk2sX78egD//+c8ArfYbGxvLzp07ee6553j88ce55pprePHFF0+6rN8VfrVwirX2b8DfAHJzc623+g0NDWN48ZeYHHWYmccX0VhfRePnfkhDUhZltWUcrj7MjqM7WHV4FQvzFpISnsJ3JnyHizIv8lYJIiI9bvHixYwbN45Ro1zBXldXx913331SG09h6EtjxowhLy+P5557jksuuaTTx5uNHj2au+++m+9973vMmTOH6dOn8+GHH3LFFVcQHh4OwJVXXsmyZcs6HNitaa3fuXPnkpmZSU5ODgATJkwgLy+vS+dqyRuBXQAMbPHa6X7PU5t8Y0wgEA2UeOHcHeIICCIkZjbLaur57i0/gH9dCot+Dv+7HKJSTrRrbGpkeeFyHlv3GPcsvYcPCz7k/in3ExTQ8aXjRETa095I2FfWr19/IqwKCwuJiIhg6tSpJ7Xx5gj7j3/8I48//jgACxYsICUlxWO7uXPn8t3vfpclS5ZQUnJ6NLR3HGDYsGGsXbuWBQsW8IMf/IBZs2YRHR3dbo2BgYE0NTWdeN3VpWNDQkJO/BwQEODVS+LeuIe9ChhqjMk0xgQDXwJeO6XNa8CN7p+/CLxnrfXaCLojMhLC2HukChueCF96Fhpq4aXboEUZAY4AZjhn8Oylz3LbmNt4dfer3LXkLuob67uzVBERnwgODqagwDWeuu+++6irqzutzbJly1i/fv1pX2dyOfyb3/zmic+3FtbgmnF+//33M3r06DM6Dq5fQMLCwrjhhhu45557WLt2LdOnT+eVV16hurqaqqoqXn75ZaZPn37S5/r3709RURElJSXU1tbyxhtvnDgWGRnp8ZeXjvTrC10ObGttA3A78DawFfivtXazMeYnxpi57mb/AOKNMbuAu4DTHv3ytYz4cMprGiitroeEofD5n0HeMtj04mltAx2B3DHuDr4/+fssyV/Cwysf7u5yRUS8bt68eSxdupSsrCzGjh3LlClTuPPOO8+or+uuu44pU6awfft2nE4n//jHP864LqfTybe+9a0zPg6wceNGJk2aRE5ODj/+8Y/5wQ9+wPjx47npppuYNGkSkydP5tZbbz3tcnhQUBA/+tGPmDRpEhdccAHDhw8/cSw+Pp6pU6eSnZ3NPffcc+L9jvTrC6abB7odlpuba1evXu21/t7dephbnlzNS984l/FpsdDUCH87D46Xwu2rISjU4+d+s+Y3PLHpCR6a/hBzBnX92T8R6Zu2bt3KiBEjeroM8SOe/k4YY9ZYa3M9te8zK51lJLgmB+QdqXK94QiACx6EsgPw6bOtfu5b475FTmIOv1j5C47VHOuOUkVERE7TZwJ7YGwYDtMisAEGzYSU8fDRH1wjbg8CHYHcP+V+Kusq+e2a33ZLrSIiIqfqM4EdHOggNbYfe0taLJJiDEy7E47ugW1vtvrZIbFDmDdiHq/ufpW9ZXu7oVoREZGT9ZnABtfEs30lVSe/OXwORDlh7ZOeP+R2c/bNhASE8JdP/+LDCkVERDzrU4GdmRDuerSr5UQ7RwDkzINd70JZfqufje8Xz7VZ17IwbyEHKw92Q7UiIiKf6VOBnREfTkVNA0erTnn2cNz1gIX1rU8+A7huuGsN2/k75vuoQhEREc/6VGBnNs8UP/WyeGwGZJ4H656GFivenColIoUZzhm8uPNF6hpPX3BARETEV/pUYDc/2rX3iIfduXLmQel+KFjTZh/XZl3L0ZqjLM1f6osSRUREPOpTge2M7UeAw5z8aFezrIshIBi2vNJmH+cMOIe40DgW7F3goypFRPzbgQMHOP/88xk5ciSjRo3id7/7XU+X1Cf0qcAOCnCQHhfG7uLK0w+GRsPgz8GWV09aX/xUgY5ALsq4iA8OfEBF3elrzIqInO0CAwP59a9/zZYtW/jkk0/44x//yJYtW3q6rLNenwpsgMFJEewq8hDYACMvc618VrC2zT4uGXQJdU11vH/gfR9UKCLiO5s3b2b27NkMGzaMBx98kDvuuINVq1Z1qo8BAwYwfvx4wLVBxogRI05sKiK+41f7YXeHIUkRLNleRENjE4EBp/y+knUxOIJcl8WdE1rtY0zCGJL6JbHkwBLmDp7bajsREY/euhcObfRun8mj4eK2Nyqqqanh6quvZv78+QwaNIjhw4czYcIEJk6ceKJNZ7fXzMvLY926dUyePLnrfwZpU98L7MQI6hst+45WMzgx4uSD/WIhYxrseBsufLDVPowxzBg4gwV7FlDXWEdwQLCPqxYR6brFixczbtw4Ro1y7cddV1fH3XfffVKbZcuWdbi/yspKrrrqKh599FGioqK8Wqucru8FdpIrpHcVVZ4e2ABDL4S374Njea7HvVox0zmTF3a8wOpDqzk39VzfFCsiZ6d2RsK+sn79+hPbQBYWFhIREcHUqVNPatPREXZ9fT1XXXUV119/PVdeeaVvCxegDwb2oETXo127iir5/CgPDZoDe+cimPTVVvuZPGAyoQGhfJD/gQJbRHqF4ODgE/ea77vvPurqTl9PoiMjbGstt9xyCyNGjOCuu+7yep3iWZ+bdBYZGkRyVCi7W5t4Fj/YNbLetbjNfkIDQxnffzwrD630fpEiIj4wb948li5dSlZWFmPHjmXKlCnceeedne5n+fLlPP3007z33nvk5OSQk5PDggV61NXX+twIG1yXxT0+2gWuHbyGXghrn4b6GggKbbWfScmTeHTtoxw5foSEfgk+qlZExDucTidr1rS9OFRHTJs27eQ9GaRb9LkRNjQHdlXrf+GGXggNx2Hfh232Myl5EgCrD632dokiIiIn6ZOBPTgpgsraBg6V13hukDENAkNd97HbMCJ+BOFB4bosLiIiPtc3A7vFxDOPgvpB+rmwZ0mb/QQ6Asntn6vAFhERn+uTgd3y0a5WZZ4Hxdug4nCbfU1Mnsi+8n0cqjrkzRJFRERO0icDOzEihKjQwHYCe4bre17bjzjkJucCsL5ovbfKExEROU2fDGxjDEPaWlMcYMBYCImGvW1vozksdhghASFsOLLBy1WKiIh8pk8GNrTzaBeAIwAyprYb2EGOIEbGj2RDsQJbRER8p88G9rD+kRyprKOksrb1Rpkz4NheKN3fZl9jEsawtWQr9Y31Xq5SRETEpc8GdlZyJADbD7Wxp3Xzfey9bd/HHp04mrqmOrYf2+6t8kRE/F5jYyPjxo1jzpw5PV1Kn9BnA3t4smtnma1tBXbiCAhLaPey+NjEsQC6LC4ifcrvfvc7RowY0dNl9Bl9NrATI0NIiAhm+6Hy1hs5HK7nsfd91GZf/cP6k9gvkY1HvLy/rYiIl23evJnZs2czbNgwHnzwQe644w5WrVrV6X7y8/N58803ufXWW31QpXjSJ9cSb5aVHMm2tkbY4Arsra9BWT5EOz02McYwJnGMRtgi0iG/WPkLth3d5tU+h8cN53uTvtdmm5qaGq6++mrmz5/PoEGDGD58OBMmTGDixIkn2nR0e80777yTX/7ylx7bim/06cAenhzFv1fso7HJEuAwnhulneP6vv8TGP3FVvsaETeCd/e/S1V9FeFB4T6oVkSkaxYvXsy4ceMYNcq1t3BdXR133333SW06sr3mG2+8QVJSEhMmTGDJkiW+KFU86NOBnZUcSU19E/tKqhiUGOG5Uf/REBwB+z9uO7DjXfdxdhzbwbikcb4oV0TOEu2NhH1l/fr1jBvn+vepsLCQiIgIpk6delKbjoywly9fzmuvvcaCBQuoqamhvLycG264gWeeecb3f4g+rE8H9gj3xLPthypaD+yAQHBOdI2w25AVmwXA1pKtCmwR8UvBwcEUFBQAcN9991FXV3dam46MsB966CEeeughAJYsWcIjjzyisO4GfXbSGcDQ/hE4TDszxcF1H/vwZjhe2mqTpLAkYkNi9WiXiPitefPmsXTpUrKyshg7dixTpkzhzjvv7OmypIP69Ag7NCiAjITwtmeKg/s+toUDK2HYhR6bGGPIisvy+kQSERFvcTqdrFmzxqt9zpw5k5kzZ3q1T/GsT4+wAYZ3ZKZ4ai44AmF/2493jYgbwc5jO6lv0opnIiLiXQrs5Cj2H62mqrah9UbBYTAgp/372HFZ1DfVs7dsr5erFBGRvq7PB3ZWciTWwo7D7Yyy086BgjVQX9Nqk+FxwwF0WVxERLyuzwd280zxdi+Lp02Bxjo42Pq+1xlRGQQ7gtlxdIc3SxSRs4S1tqdLED9xJn8X+nxgD4zrR2RoIJsKytpu6HSvBJTf+hJ+AY4AMqMz2V2224sVisjZIDQ0lJKSEoW2YK2lpKSE0NDQTn2uT88SB9fs7uyUaDYVtjNTPLI/xKS1GdgAg2MGs65onRcrFJGzgdPpJD8/n+Li4p4uRfxAaGgoTqfn5a5b0+cDGyA7NYonP95HfWMTQQFtXHTowAIqg2MGs2DvAi1RKiInCQoKIjMzs6fLkF6sz18SB8hOjaauoYldRZVtN3ROgvICKCtotcngmMEA7Cnd480SRUSkj1Ng4wpsgI1euI89JGYIgO5ji4iIVymwgcz4cMKDA9jcXmAnj4aAkDYD2xnhJNgRzO5SBbaIiHiPAhtwOAyjUqLbH2EHBkNKTsdmiiuwRUTEixTYbqNSo9hysJzGpnYeuXBOhML10HD6LjfNBsUMUmCLiIhXKbDdRqdGU1PfxO7i9iaeTYTGWji8sdUmQ2KGUFhVSHV9tZerFBGRvkqB7dY88azDC6gcaP2y+ImZ4mWaKS4iIt6hwHYbnBhBaJCj/fvY0akQmdLmfezB0a7A3lW6y5sliohIH6bAdgtwGEYOiGp/hA0wcGKbgT0wciBBjiA9iy0iIl7TpcA2xsQZYxYZY3a6v8d6aJNjjPnYGLPZGLPBGHNtV87pS2OcMWwqKKehsanths6JULoPKos8Hg5wBJAWmca+8n0+qFJERPqiro6w7wXetdYOBd51vz5VNfBla+0o4CLgUWNMTBfP6xPj0mI4Xt/I9va22nROcn0/sLLVJulR6QpsERHxmq4G9mXAk+6fnwQuP7WBtXaHtXan++dCoAhI7OJ5fWLcQNcFgvUHSttuOGAsOAKhYHWrTdKj09lfsZ/GpkZvligiIn1UVwO7v7X2oPvnQ0D/thobYyYBwYBfPqQ8MK4fceHBrNvfTmAHhbpWPctvPbAzozKpb6qnsKrQy1WKiEhf1O5uXcaYxUCyh0Pfb/nCWmuNMa2uOmKMGQA8DdxorfV4k9gYcxtwG0BaWlp7pXmdMYZxA2NYt/9Y+41Tc+HT56CpERwBpx1Oj0oHYF/5PgZGDvR2qSIi0se0O8K21s621mZ7+HoVOOwO4uZA9jgLyxgTBbwJfN9a2+r+lNbav1lrc621uYmJPXPVPGdgDLuLqyg7Xt92Q2cu1FVC8TaPh1sGtoiISFd19ZL4a8CN7p9vBF49tYExJhh4GXjKWvtCF8/nc+PSXPexP23vPnZqrut7K5fF40LjiAyKJK8sz4vViYhIX9XVwH4YuMAYsxOY7X6NMSbXGPN3d5trgBnATcaY9e6vnC6e12fGDIzGmA5MPIsfDKExrU48M8ZopriIiHhNu/ew22KtLQFmeXh/NXCr++dngGe6cp7uFBUaxJDEiPbvYxvjuiyev6bVJunR6aw7vM7LFYqISF+klc48GJcWw/oDpVjbzs5dqblQvBVqPT+3nR6VzsGqg9Q01PigShER6UsU2B7kDIzlWHU9+0ra2W3LmQu2CQo9j6IzojKwWA5UHPBBlSIi0pcosD0Yn+5aiG31vnYui6dOcH1vZeKZZoqLiIi3KLA9GJYUSVRoIKv2Hm27YVgcxA2GAs/3sZsDO688z8sViohIX6PA9sDhMORmxLEqr53ABvfEs1Xg4X53eFA4Cf0S2F++3wdViohIX6LAbsXEjDj2HKmiuKK27YapuVB5GMryPR4eGDlQ97BFRKTLFNitmJQZB8Dq9kbZTvcCKq08j63AFhERb1Bgt2J0ajShQQ5WthfY/bMhIKTViWfOSCdF1UXUNrYzUhcREWmDArsVwYEOcgbGtH8fOzDYtd1mKxPPnBFOLJaCigIfVCkiIn2FArsNkzLi2FJYTkVNBzYCKVwHjae3a96pK7/S8z1uERGRjlBgt2FiZhxNFta2tz926gRoqIHDm0871BzYuo8tIiJdocBuw/i0WAIcpv3nsZ0TXd89TDyLC40jLDBMgS0iIl2iwG5DeEgg2SlRrNhb0nbDmDQIT/S4EYgxRjPFRUSkyxTY7ThncDzrD5RSXdfQeiNjXM9j56/yeFiBLSIiXaXAbsfUwQnUN1pWtntZfAKU7ITjp68/PjByIAUVBTTZJh9VKSIiZzsFdjsmZsQRFGD4aHc7l8VP3Mdee/qhSCd1TXUUVRf5oEIREekLFNjt6BccwLi0WD7afaTthinjAePxeWzNFBcRka5SYHfA1MEJbC4sp7S6rvVGoVGQmOVxxTMFtoiIdJUCuwOmDonHWvi4vcviqZ537koOTybQBJJfocVTRETkzCiwO2CMM4aw4IAO3MfOheNH4djek94OdAQyIGKARtgiInLGFNgdEBzoYFJmHMvbu4/dvHOXh+ex9WiXiIh0hQK7g6YOTmBPcRUHy4633ihxBASFeVzxTIEtIiJdocDuoOnDEgBYuqO49UYBgZAyzuMCKgMjB1JeV05ZbZmvShQRkbOYAruDsvpHkhwVypLtbQQ2uC6LH9oIDSfvf+2MdAJo4pmIiJwRBXYHGWM4b1giH+48Qn1jGyuWpeZCY50rtFvQo10iItIVCuxOmJmVSEVtA+va2m7zxMSzk+9jOyPcI2ztiy0iImdAgd0JU4cmEOAwfLCjjSVGo1IgMuW0+9hhQWHEh8brkriIiJwRBXYnRIUGMSEttmP3sT3MFHdGOjXCFhGRM6LA7qTzshLZXFhOUUVN642cuXAsD6pOfm7bGenUCFtERM6IAruTZmYlArB0RxuLqKS672OfshGIM8LJoapD1DfV+6o8ERE5SymwO2nkgCiSIkN4b9vh1hul5IAJOO0+tjPSSaNt5FDVIR9XKSIiZxsFdicZY5g1oj8fbC+mtqHRc6PgcOg/8rSZ4qkRqYCexRYRkc5TYJ+BC0f1p6quse3NQFJzoWAtNH32zHbzs9iaeCYiIp2lwD4D5w6OJzw4gHc2t3FZ3JkLtWVQsuvEW4n9EglyBGmELSIinabAPgMhgQHMzEpi8dbDNDVZz41OTDz77LJ4gCOA1IhUBbaIiHSaAvsMXTCyP8UVtXya38qqZwnDICTqtIlnqZGpuiQuIiKdpsA+Q+dnJRHgMLyzpZXL4g4HpI73uESpRtgiItJZCuwzFB0WxDmD4ljUWmCD67L44c1QV33iLW2zKSIiZ0KB3QUXjOjPrqJK9hRXem7gzAXbCAc//ewt9yYgBZUF3VGiiIicJRTYXXDBqGQAFm5uZSGU5olnLe5ja19sERE5EwrsLkiN6UfOwBje3HDQc4OIRIhJP2mm+InFUzTxTEREOkGB3UVzxgxgc2E5e49UeW7gzIX8z9YUjwiOIDYkViNsERHpFAV2F10yegAAb3xa6LlBai6U50P5Z6Nw7dolIiKdpcDuopSYfuSmx/JGa5fFnacvoOKM0L7YIiLSOQpsL5gzZgDbD1ew83DF6QeTx4Aj6KTnsZ2RTg5WHqShqaEbqxQRkd5Mge0Fl4wegDF4HmUHhULy6JP2xnZGOmmwDRyubuMZbhERkRYU2F6QFBXK5Mw43thQiLUe1hZ3Nu/c5dqOs/lZbN3HFhGRjlJge8mlY1LYXVzF1oMeLoun5kJ9FRRtBfQstoiIdJ4C20suHT2AQIfh5XUeQviUiWf9w/oTaAI18UxERDpMge0lcdy9ixMAACAASURBVOHBnD88iZfXFdLQ2HTKwUHQL+7ExLMARwApESkaYYuISIcpsL3oqvFOjlTWsmzXkZMPGAOpE06beKbAFhGRjupSYBtj4owxi4wxO93fY9toG2WMyTfGPNaVc/qzzw1PIjYsiBfXeLosPtF1D7vGtUuXnsUWEZHO6OoI+17gXWvtUOBd9+vWPAgs7eL5/FpwoIO5Y1N4Z8thyo7Xn3ww7RzAwoGVgGuEXVpbSmVdKzt9iYiItNDVwL4MeNL985PA5Z4aGWMmAP2Bd7p4Pr935XgndQ1Np28I4pwIjkDY95HrZaS22RQRkY7ramD3t9Y2J9MhXKF8EmOMA/g18N32OjPG3GaMWW2MWV1cXNzF0nrGGGc0Q5IieHHtKZe7g8NgQA7s/xhosWuX7mOLiEgHtBvYxpjFxphNHr4ua9nOulYM8bBqCN8AFlhr200ma+3frLW51trcxMTEDv8h/IkxhqvGO1mz7xi7i0+53J0+xTXxrL7ms2exdR9bREQ6oN3AttbOttZme/h6FThsjBkA4P5e5KGLKcDtxpg84BHgy8aYh734Z/A7V01IJdBh+M/K/ScfSDsXGuugcC1RwVFEBUdxoOJAzxQpIiK9Slcvib8G3Oj++Ubg1VMbWGuvt9amWWszcF0Wf8pa29bktF4vKTKUC0f154U1+dTUN352IO0c1/cW97E1whYRkY7oamA/DFxgjNkJzHa/xhiTa4z5e1eL683mTUrnWHU9Czcd+uzNsDhIHHHiPrYzwklBhSadiYhI+7oU2NbaEmvtLGvtUPel86Pu91dba2/10P5f1trbu3LO3uLcwfGkx4fx7IpTLounT3E92tXUiDPSSUFlAY1NjZ47ERERcdNKZz7icBium5TGyryjJ++TnXYu1JbD4U04I53UN9VTfLx3zogXEZHuo8D2oS9OcBIUYHi25eSz9Cmu7/s+PrHNpiaeiYhIexTYPpQQEcLnRyXz4pp8qusaXG9GOyE6DfZ/pG02RUSkwxTYPnbjuRmU1zTw0toWk8vSp8C+j0kO60+ACdBMcRERaZcC28dy02MZnRrNE8v30tTkXlcmbQpUFRFUeoDk8GSNsEVEpF0KbB8zxnDztAx2F1d9tu1m+rmu7/s+0rPYIiLSIQrsbnDp6BQSI0P454d7XW8kDIOweFdgR2hfbBERaZ8CuxsEBzr4n3PS+WBHMbuKKsAYSJ8KectwRjg5WnOU6vrqni5TRET8mAK7m8ybnEZwoIMnlue53sicAWUHcDpCAG0CIiIibVNgd5OEiBCuyEnlhTX5FFfUugIbGFhaCOjRLhERaZsCuxvddt4g6hqbeGL5Xtd97Ij+OA9tBRTYIiLSNgV2NxqcGMHF2ck8/fE+ymsbIGM6UXkfERkUqUviIiLSJgV2N/vGzCFU1DbwzCf7IHMGpvIwzn4JGmGLiEibFNjdLDs1mhnDEvnnh3upHTgVAKd1aIQtIiJtUmD3gG/MHMyRyjqe3xUA0QNxVldQUFFAk23q6dJERMRPKbB7wOTMOCakx/LXpXtpSp+G8+h+6prqKK7WNpsiIuKZArsHGGO4/XNDKCg9zgqycR4vB/QstoiItE6B3UNmDktkfFoMD29LxFnv2npTE89ERKQ1CuweYozhuxdm8Wl5BOFBAzBohC0iIq1TYPegc4ckMGVQPB/UDCe5oZH88v09XZKIiPgpBXYPu/vCYSyqHYmzvp6Cku09XY6IiPgpBXYPy82Ig0HnkdrQwIHKgp4uR0RE/JQC2w/870UTcNRFc6SphuMNx3u6HBER8UMKbD8wxhlDSL8sADbsXdvD1YiIiD9SYPuJqRMuBWDBB/N7uBIREfFHCmw/MTr7YgCqjq1jY35ZD1cjIiL+RoHtJ2L7JRCGg/CQQ/z0jc1Ya3u6JBER8SMKbD9hjMEZGs/RwAaO7dvAO1sO93RJIiLiRxTYfsQZO4T8oECuitrGT9/cQk19Y0+XJCIifkKB7UecsUPJDwrmmrgdHDh6nD+9v6unSxIRET+hwPYjzkgnNQYaj6zli2Pi+MsHe8g7UtXTZYmIiB9QYPsRZ4QTgHxHIz8YeYTgQAf3v6YJaCIiosD2K85IV2AfCA0nJv89vnPBMD7YUczbmw/1cGUiItLTFNh+JCUiBYMhP3EI7HibG89JY3hyJA+8toXymvqeLk9ERHqQAtuPhASEkBSWRH5UIpQXEFi8mYeuHE1RRQ0PLdjW0+WJiEgPUmD7GWekk/wAAxjYsZBxabHcMi2T51bu56PdR3q6PBER6SEKbD/jjHCSX30YnLmw/S0A7rogi4z4MO59cSPVdQ09XKGIiPQEBbafcUY6KaouonbIbChcCxWH6RccwMNXjWH/0Wp+/c6Oni5RRER6gALbzzTPFM935rje2Pk2AOcMiueGc9L45/K9rNl3tKfKExGRHqLA9jMZURkA5AUFQfRA2L7wxLF7Lx5BSnQ/vvP8p1TW6tK4iEhfosD2M+lR6QDklefBsM/Dnveh/jgAESGB/PbaHPKPVfPj1zb3YJUiItLdFNh+JjI4kvjQeFdgD58D9dWw690TxydlxvGNmUOYvyafBRsP9lyhIiLSrRTYfigjOoN95fsgYxr0i4Wtr590/NuzhzLWGc19L23kYNnxHqpSRES6kwLbD2VEZZBXlgcBQZB1ievxroa6E8eDAhw8+qVx1DU0cfd/P6WpSWuNi4ic7RTYfigzOpNjtccoqy2DEXOhtgz2Lj25TUI4D8wdyUe7S/jDe9qGU0TkbKfA9kPNM8X3lu2FQTMhOBK2vnpau2tyB3LluFQefXcHy3YWd2uNIiLSvRTYfuikmeJBoTDsQtj2JjQ1ntTOGMNPr8hmaFIE3/7Pet3PFhE5iymw/VBqZCqBJtB1Hxtcl8WrS2DfR6e1DQsO5M83TKC2vpFv/nst9Y1N3VusiIh0CwW2HwpyBOGMdLpmigMMvQAC+8HW1zy2H5wYwcNXjWHt/lJ+vmBrN1YqIiLdRYHtpzKiM1yXxAGCw2HILNjy2mmXxZt9YWwKX5mawRPL8/jvqgPdV6iIiHSLLgW2MSbOGLPIGLPT/T22lXZpxph3jDFbjTFbjDEZXTlvX5ARlcH+8v00Ngd09lVQeQj2LW/1M9+/ZATThiTw/Vc2sjpP642LiJxNujrCvhd411o7FHjX/dqTp4BfWWtHAJOAoi6e96yXEZVBXVMdhVWFrjeGXQTBEbBxfqufCQxw8Md543HGhvH1Z9ZQUKpJaCIiZ4uuBvZlwJPun58ELj+1gTFmJBBorV0EYK2ttNZWd/G8Z72M6AyAzyaeBYfB8Etdl8Ubalv9XHRYEI9/OZfahiZufXK19s8WETlLdDWw+1trmxe0PgT099BmGFBqjHnJGLPOGPMrY0yAp86MMbcZY1YbY1YXF/ft54pP7NrVfB8bIPuLUFN60tringxJiuAP141j+6Fy7nh2HQ2aOS4i0uu1G9jGmMXGmE0evi5r2c5aawFPa2QGAtOB7wITgUHATZ7OZa39m7U211qbm5iY2Nk/y1klLjSOyOBI1+IpzQafD/3iYNML7X5+ZlYSP7ksm3e3FfHDVzfh+p9HRER6q8D2GlhrZ7d2zBhz2BgzwFp70BgzAM/3pvOB9dbaPe7PvAKcA/zjDGvuE4wxDI4ezO7S3Z+9GRAEoy6H9c9BbSWERLTZxw3npHOorIbH3t9FclQ/vj17qI+rFhERX+nqJfHXgBvdP98InL5+JqwCYowxzUPmzwFbunjePmFwzGB2l+0+eXQ8+mpoOO7aEKQD7r5wGFeNd/LbxTv4z8r9PqpURER8rauB/TBwgTFmJzDb/RpjTK4x5u8A1tpGXJfD3zXGbAQM8HgXz9snDIkZQlltGSU1JZ+9OfAciHLChuc71IcxhoevGs2MYYl8/5VNLNpy2EfVioiIL3UpsK21JdbaWdbaodba2dbao+73V1trb23RbpG1doy1drS19iZrbV3rvUqzwTGDAU6+LO5wwNhrYfe7UH6wlU+eLCjAwZ+vH092ShTf/PdabRQiItILaaUzP9Yc2LtKT9k+M+d6sE3w6XMd7is8JJAnb57EoMRwvvrUalbsKWn/QyIi4jcU2H4ssV8ikcGRJ4+wAeIHQ9oUWP9v6MTs75iwYJ65dTLO2DBu/tcq1u0/5uWKRUTEVxTYfswYw5CYIacHNrhG2SW74MDKTvWZEBHCv2+dTEJkCDf+cyWbCsq8VK2IiPiSAtvPDY4ZzK7SXac/Rz3qcggKg/XPdLrP/lGh/PvWyUSGBjHv8U9Yf6DUS9WKiIivKLD93JCYIZTXlZ88UxwgJBJGXg6bXoa6qk7364wN4/mvnUNseDA3/H0Fq7RZiIiIX1Ng+7lWJ54BjLse6ipg6+tn1LczNoznb5tC/6gQvvyPlSzfdaQrpYqIiA8psP3ckJghAJ7vY6dPhbjBsPqJM+4/OTqU/9w2hfT4ML7yr1W8t03PaYuI+CMFtp+LD40nOiTa8wjbGMj9Chz4BA5tOuNzJEaG8NxXzyGrfyRffWoN81cf6ELFIiLiCwpsP9e8pvie0j2eG+RcD4GhsLprS7PHhgfz3G3ncO7geO55YQOPvbdTG4aIiPgRBXYv0OpMcYCwOBh1JWz4L9RWdOk8ESGB/OPGiVwxLpVH3tnBj17dTGOTQltExB8osHuBobFDKa8rp6ja02ZowMRboK6yw+uLtyU40MGvrx7L184bxNOf7OMb/17D8brGLvcrIiJdo8DuBbJiswDYfmy75wapE2DAWFj1j06tfNYah8Nw38UjuP8LI3lny2Gu+evHHCqr6XK/IiJy5hTYvcCw2GEAbD/aSmAbA7m3QNEW2P+x1877lamZ/P3LuewprmTuYx9qgRURkR6kwO4FIoIjcEY42XZ0W+uNRn8RQmPgkz959dyzRvTnpW9MJTjQwbV//ZhX1xd4tX8REekYBXYvkRWXxY5jO1pvEBwOuTfD1jfgaCszys/03MmRvPrNqYx1xvDt/6zn4be20dDY5NVziIhI2xTYvURWXBb7yvdRXV/deqNJt4EjED75i9fPHx8RwjO3Tmbe5DT+8sFubvjHCooqdF9bRKS7KLB7iazYLCyWnaU7W28UNcB1aXzdM3Dc+1tnBgc6+PkVo/n11WNZf6CUOb//kJV7tQa5iEh3UGD3EsPjhgNtTDxrNuWbUF8Fa/7ls1qumuDklW9OJTwkkOse/4S/Ld2tRVZERHxMgd1LDAgfQGRwZPuBnTwaBs2EFX+Fhjqf1TM8OYrXbp/KhSP78/MF27jlydUcqaz12flERPo6BXYvYYwhKzaLbcfamCnebModUHEQNv7XpzVFhgbxp+vH88AXRvLhriNc9OgylmxvZXEXERHpEgV2L5IVl8XOYztpbGpn5bEhsyB5DCz7NbTXtouMMdw0NZPXb59GfHgwNz2xigde20xNvVZHExHxJgV2LzI8bjjHG46zr3xf2w2NgRn3uB7v2vxyt9SWlRzJq7dP5aZzM/jXR3lc9thyNheWdcu5RUT6AgV2L5Idnw3AppIObKU5fA4kDoelj0BT9zwzHRoUwANzR/HEVyZytLqOyx5bzm/e2U5dg57ZFhHpKgV2L5IZnUm/wH5sPrK5/cYOB0z/LhRvhW1v+L64Fs7PSmLRd2Ywd2wKv39vF1/4w4dsyNeypiIiXaHA7kUCHAGMiBvRsRE2wKgrIG4QLP2VVzYF6YyYsGB+c20O/7wpl9LjdVz+x+X8YuE23dsWETlDCuxeJjshm+1Ht1PfVN9+44BA1yj70AbY+rrvi/Pgc8P78853zuOLE5z8ecluLv7dMpbtLO6RWkREejMFdi+TnZBNbWMtu0t3d+wDY66FhGHw3k99PmO8NdH9gvjlF8fy9C2TsNbyP/9Yye3PruVwuZY2FRHpKAV2LzMqfhQAm4508LJ4QCCc/304sh02PO/Dyto3fWgiC++cwXdmD+OdLYeZ9esP+OeHe7WRiIhIByiwe5mBkQOJCo7qeGADjLwMBuTA+w9BQ8+uRhYaFMC3Zw/lnTtnMD49lp+8sYUvPLacj3Yf6dG6RET8nQK7lzHGMCp+FFtKtnTmQzDrR1C236drjHdGRkI4T35lIn+6fjzlx+uZ9/gKvvrUavYeqerp0kRE/JICuxfKTshm57Gd1DR04h7w4M9BxnT44JdQ4x8LmhhjuGT0AN69+zzu+XwWH+06wgW/+YCfvL6FsuoOTKoTEelDFNi90OiE0TTYhs6Psi98EKpLXEuW+pHQoAC+ef4Q3r9nJl+c4OSJj/Zy3iPv8/dle/QYmIiImwK7F8pJygFgXdG6zn0wZRzkzINP/uxattTPJEWG8vBVY3jzjulkp0Tz0ze3cv4jS3hu5X7qNTFNRPo4BXYvFBsaS0ZUBuuL1nf+w7N+BI4geOeH3i/MS0amRPHMrZN59tbJJEeHct9LG7nwt0t57dNCmpq077aI9E0K7F4qJymH9cXrsZ1dwSwyGabf5VqudO9S3xTnJecOSeCl/z2Xv385l5BAB996bh2X/H4Zb28+pOAWkT5Hgd1LjUsaR2ltKXnleZ3/8JTbISYN3voeNPr35C5jDLNH9mfBt6bz++vGUVPfyNeeXsPFv1vGq+sLaFRwi0gfocDupZrvY5/RZfGgULj4l1C0BT5+zMuV+YbDYZg7NoXFd53Ho9fm0Ggt3/7Pemb/5gP+u+qAdgQTkbOeAruXyojKIDokuvMTz5plXQwjvgBLfgFH93q3OB8KDHBw+bhU3rlzBn+5YTxhwQH834sbOP+RJTz1cR7H6zSrXETOTgrsXsphHOQk5px5YINrlO0IhDfv7vbdvLrK4TBclD2AN+6YxhNfmUhydCg/enUzUx5+l1+9vU3rlIvIWUeB3YvlJOWQV57H0ZqjZ9ZBVArM+iHsfhc2v+Td4rqJMYbzs5J44etT+O/XpjA5M44/LdnNtF+8x13/Xc/mQv9YJEZEpKsU2L1Ybv9cAFYfWn3mnUy81fV89lv3QvUZBr8fMMYwKTOOv/5PLku+O5PrJ6ezcNMhLv39h8x7/BMWbTmsCWoi0qspsHuxUQmjCAsMY+WhlWfeiSMA5v4Bjh+DN+/yXnE9KD0+nAfmjuLje2dx38XD2Xukiq8+tZrpv3iPP7y7k6IKXS4Xkd5Hgd2LBTmCmNB/QtcCGyB5NMy8Fza/DBtf8E5xfiA6LIivnTeYZf93Pn+5YQKDEiP49aIdnPvQe3zz2bV8vLuk88+xi4j0kMCeLkC6ZvKAySxbvYyi6iKSwpLOvKOpd8KOha4JaOlTIWqA94rsYYEBDi7KTuai7GT2FFfy7Ir9zF+Tz5sbDjIkKYIvTRzIFeNSiY8I6elSRURapRF2LzcpeRJA10fZAYFwxV9d+2W/dnuvmzXeUYMSI/jBnJGs+H+z+NUXxxAREshP39zK5J+/y21PrWbRlsNat1xE/JJG2L1cVlwWUcFRrDy4kjmD5nSts/jBcMFP4K17YMVf4Jz/9U6Rfig0KICrcwdyde5Adhyu4IU1+by0toB3thwmISKYK8alcnXuQIb1j+zpUkVEADD+eg8vNzfXrl7dhdnPfch33v8OW49uZeFVC7vembXw3HWwazHc8jakTuh6n71EfWMTH2wvZv6aA7y7tYiGJsvo1Gguy0lhzpgUkqNDe7pEETnLGWPWWGtzPR5TYPd+z217jp+v+DkLrljAwKiBXe+w+ij8dYZrD+2vLYN+MV3vs5cpqazllfWFvLKugI0FZRgDkzLimJuTwiXZA4gND+7pEkXkLNRWYOse9llgaspUAJYVLPNOh2Fx8MV/QnkhvPrNs/Z+dlviI0K4ZVomr98xjffuPo87Zw2juLKW77+8iYk/W8zN/1rFK+sKqKxt6OlSRaSP0Aj7LDHn5TkMjBzIn2f/2XudfvQHeOcHcOHP4NzbvddvL2WtZXNhOa9/WsjrnxZSWFZDcKCDGUMTuHBUMheM6K+Rt4h0SVsj7C5NOjPGxAHPAxlAHnCNtfaYh3a/BC7FNaJfBHzb+utvCr3U9NTpzN8xn+MNx+kX2M87nU65HfZ/Aot+CEkjYMgs7/TbSxljyE6NJjs1mu9dNJw1+4+xcNMhFm46xOKtRQQ4DJMz47goO5kLRybrnreIeFWXRtjuID5qrX3YGHMvEGut/d4pbc4FfgXMcL/1IXCftXZJW31rhN05HxV+xNcWfY0/zvojM5wz2v9AR9VWwj8uhPJ8+Or7rpnkcpLmkffCTYd4a9NBdhdXATAuLYbPj0pm9ogkBidGYIzp4UpFxN/5bNKZMWY7MNNae9AYMwBYYq3NOqXNFOAxYBpggKXA/1hrt7bVtwK7c+oa65j2n2lcNvgyvn/O973b+bF98Pj5EBYPty6G0Gjv9n+W2VVU4Rp5bz7EpoJyAAbG9eNzWUmcPzyJcwbFExoU0MNViog/8mVgl1prY9w/G+BY8+tT2j0C3IorsB+z1rabKArszrvjvTvYeWwnb135lvdHc3kfwlOXQeZ5MO95CAjybv9nqcLS47y/vYj3txXx4a4j1NQ30S8ogKlD4pmZlcTnhieREuOlWxgi0ut16R62MWYxkOzh0Emha621xpjT0t8YMwQYATjdby0yxky31p42pdkYcxtwG0BaWlp7pckppqdOZ8mBJewu3c2Q2CHe7TxjGsz5Lbx2B7z2Lbj8T67HvqRNKTH9uH5yOtdPTqemvpGP95Tw/rYi3ttWxOKtRQAMT45k2pAEpg1NYHJmPP2CNfoWkdN1xyXxe4BQa+2D7tc/Amqstb9sq2+NsDuvqLqI2fNn842cb/D1sV/3zUmWPAxLHoLpd8OsH/nmHH2AtZZdRZW8t62ID3YUszrvGHWNTQQHOBifHsP0oYlMG5JAdmo0AQ79YiTSV/jykvivgJIWk87irLX/d0qba4GvAhfhuiS+EHjUWvt6W30rsM/Ml9/6MlX1Vbw490XfnMBaeONOWPMvuPhXMPk235ynjzle18jKvKMs33WEZTuPsPWg6953dL8gzh0cz7ShCUwZFE9mQrgmr4mcxXz2WBfwMPBfY8wtwD7gGvcJc4GvW2tvBV4APgdsBCywsL2wljN3QfoF/HLVL9lXvo/0qHTvn8AYuOTXUFkEb/0fhERAzjzvn6eP6RccwHnDEjlvWCIAxRW1fLT7CB/uPMKHu47w1qZDACRGhjApM45zMuOYPCieoUmafS7SV2jhlLPMwcqDXPjihXx7/Le5dfStvjtR/XF47kuw5wO48m8w5hrfnauPs9ay50gVK/YcZcXeElbsOcqh8hoA4sKDmZQRx+RBcUzKjGNEchQOXUIX6bW0lngfM+/NeTTaRp6f87xvT1RXDc9eA/uWw1V/h+yrfHs+AVwBvv9oNSv2Hj0R4vnHjgMQFRrIhPRYxqfFMj49lrEDY4gI0aZ8Ir2FLy+Jix+6IP0CfrPmN+RX5OOMdLb/gTMVHOZ6xOvfV8OLXwXjgFFX+O58ArhWXEuPDyc9Ppxrcl2bvRSUHmfFnhJW7j3Kmn3HeH97MQAOA1nJUYxPizkR4hnxYbqMLtILaYR9FiqsLOTzL36e23Nu52tjv+b7E9ZWuEL7wAqY8yhMuNH355Q2lVXXs+7AMdbuL2Xd/mOs319KhXujkrjwYMYNjGF8eiw5A2PITo0mup+eqxfxB7ok3gfd/PbNFFUX8frlr3fPaKquCv77Zdc+2rN/DNPu9P05pcMam1yPka3Zd4y1+11fe9xLqAJkxIcx2hnDmNRoxjijGZUarUvpIj1Al8T7oLmD5/LD5T/k0+JPyUnK8f0Jg8PhS8/BK1+HxffD8aOu4NalV78Q4DBkJUeSlRzJvMmuRYlKq+vYWFDGhvwyNuSXsibvKK9/Wgi4/mcbnBjBmNRoRjtdIT5yQLQWdRHpQQrss9QF6Rfw8xU/57Xdr3VPYAMEBsOVj0NoDCz/HZQfhLl/gCDtWuWPYsKCmT40kelDE0+8V1xRyyZ3iG8sKGXZriO8tK4AcN0Pz0gIZ+SAKEYMiGLkgChGpkSRFBmie+Ii3UCBfZYKDwpnVtosFuYt5HuTvkdIQEj3nNgRAJf+GqJT4d2fwLE8+NKzEJHY7kel5yVGhnD+cNcmJc0Ol9fw6YFSNheWs/VgOesPlPLGhoMnjseFB7tDPJKRKa4wH5wYQVCAoyf+CCJnLd3DPot9XPgxty26jYenP8ylgy7t/gI2vwIvfx3CE12zyfuP7P4axCfKa+rZdrCCLYVlbD1YwZaD5Ww/XEFdQxMAwQEOhiRFMKx/BEP7RzKsfyTD+kcwMDZMz4mLtEGTzvqoJtvEF17+AnGhcTx9ydM9U0TBWnjuOqirdF0ez76yZ+oQn2tobGLvkSq2HCxny8Fyth2sYOfhCgrLak60CQ1yB3lSpDvIIxjWP5LUmH4KchEU2H3aU5uf4lerf8X8L8xneNzwnimirADm3wT5K2Hy/8IFP3Hd75Y+oaKmnp1Flew8XMGOw5XsOFzBzsOVJ1ZrA+gXFMCQpAiG9o9gSFIEgxLCGZQYQXp8GCGBmugmfYcCuw8rqy1j9vzZXDroUh4494GeK6ShDhb9CFb8GZyT4Op/ue5zS59VdryeXUWuEN95uJKdRRXsOFzB4fLaE20cBpyxYQxKDGdQQoT7uyvM+0dpspucfRTYfdwDHz3Am3veZPHVi4kOie7ZYja95NpT2xHo2l9bl8jlFBU19ew9UsXeI1XsLq5iT3Ele4pdr4/XN55oFx4cQGaLIM9MCCctLoz0+HBiw4IU5tIrKbD7uG1Ht3H161f7fkOQjirZDS/dBgWrYcy1cMmvILSHf5EQv9fUZDlUXsOe4ir2HKl0f3cFekHpcVr+UxYZEkhafBjp8WGkxYWTHh9GelwYafFhDIjupz3GxW8psIWvL/o6W49u5e2r3iY00A+ei25sgGWPwAe/hKgUsICnmAAAFVlJREFUuPxPkDmjp6uSXqqmvpEDR6vZV1LNvqPV7CupYl9JNfuPVpN/rJr/396dBsdR3nkc/z5zakajYyRZkmXZxgYbg2NibhIIdqiEQKoIAW+CoUKxFQjZqs1F4uxmIVWBJLtQWxsSXqSyIUAgtWGJCSZQhCPZcAQTMLaJDdjGXMa2ZMmydY9mRnM9++LpGfVII8lGR89o/h+qq7uf7p55/FSj33T3093J9MjfOZ/bRWtdgMXW0bg5Kg/SGg6yIByQJ7wJR0lgC7Z2buXLz3yZW869hfUr1jtdnRFt22DTV6DnfTj9Wrj4RxAIO10rMYekM5pDfTEO5AJ9iAPd0VygR6xnrGfVBr0sqA3QGg6woDZoxmEz3xoOynPXxYySwBZorbn2qWs5GjvKE1c8gcdVREcRiSg8fzu8/HMI1sOld8DKK+WxpmLGaa3pGUqwvydKe2+Mtt4Y7X1RM7bm7dfNwZxutwd4NtxbagPMr6mgIeSXW9TEhyaBLQB4/uDzfP3Zr/PDj/+QK5YV4WswO3bC49+Ajh1w0qfgM7fDvOVO10qUsWygt/fFbCEezc239cbGHKF7XIqm6gpaaitorjEhPjJIqIuJSWALwPzxueaP13A0bo6yZ+1xpccjnYJXfwnP3wHJKJxzI6z5FzlNLoqS1pqBWIqDvVE6+uN09Mfo6I/T2R/nUF+MzoE4Hf3x3BPgsrKhPr+mguaaClpqAzRb8001FTRVVzAv5Mfnkce7lhsJbJGzpWMLN/zpBjactYHrVhbxe6sjR+C5H8P2B0xYX3QLnHEduOX6oSgt2aN0E+hxOq1QzwZ8pzU9PCrUAcJBrwnvKj9N1RU0VvlpzE5X+2msMssqvPJwmblCAlvk+eqfv8qu7l08deVTVPmqnK7OxDpeh6f/DfZvhrqlsPZm+Mg6cMmRh5g7tNb0RpMc6ovRNRina2CYwwPDZnpwmK4BMz4yOEwqM/Zvdk3AS5MV4Nkgb6zy01jtpyFkhnkhP9UBj9yfXuQksEWePd17uOqJq7j21Gv57tnfdbo6k9Ma3n4anv0xHH4TGleaI+6TPysd00RZyWQ0PdGECfTBOEesUM8PdzNtv5Uty+tW1Ff6aajymXHIT0PIZ8b2siofdUEfHnnj2qyTwBZj3PbybTz6zqNsvGwjy8Ml0rErk4Fdm+C5/4Ce92D+R+GCm+CUz5nXegohgJEj9q7BON2RBEcj5ui8eyjB0cFhjkbs0wkS6bGn45WCcNBHQ8gK8io/9ZU+5lX5qav0EQ76qKv0UVfpJRz0URv0yQNppoEEthijL97HZX+4jKU1S7n/kvtL6zRZOgWvPwSbfwrd70LdiXD+N+Gj68FThB3phChiWmsGh1O58O6OmEA/YpseKU+M6RWfpRTUBryEK83ReXZcF7LNW+FeX+knXOkl5JdT9KNJYIuCHnn7EW59+dbivc1rMpk0vPUEvHinuRUs1AznfMV0TgvNc7p2QsxJ8WSa3miC7kiC3miCnqEEvUMJeqJJM7aG3LJoouDpeTCn6LNH6tlxbdBrhoCPmqCX2oCX2qCPmoAprwl453QnOwlsUVBGZ7j+mevZ07OHTZ/bREuoxekqfThaw74XzBH3+8+D22cevHLOjdB6ptO1E6Ksaa2JDKfoHUrSPTRsBbkV7lET9t250E/QH03SF0uSLtC5LqvC66I24MsFeDbga4NeK+TzA978CPBR6XMX/RG9BLYYV9tgG+seX8eqhlXcffHduFSJdzI5she23gM7HoREBFpOh7NvgFM/D/6Q07UTQhyDbMj3RZP0x5L0RZP0xRK5aTNOWOVJK+QT9EaTY+55t/O4FLVBL9UBL9UV2bHHNu+hZoJls/FudglsMaHfv/17bnv5tuK/N/t4xAfg9d/Bq7+Co3vBWwmnXg6rr4HF58ttYULMUfFkeiTko4m8QM8G/EDMhP5APMVgLMlA3MyPd+o+y+9x5YK8JmBC/esXncSZi+umrf4TBXYRPVBaOGXdsnW81P4SP93+U1Y1rOKMpjOcrtLUVVSb69ln3wAHX4UdvzXv4t75INQuNsG96gtQf6LTNRVCTKMKr5sKr5um6uN7K6HWmuFUxgS5FeIDsZQ1Hgl4+7KeoQSpSUJ+OskRtgBgMDHI+ifWE0vF2HjZRhoCDU5XafoloqaT2t//B/b9FdDQvApWXmFOmUt4CyEcJqfExTHZ27OXLz35JU6qPYl7P3MvQW/Q6SrNnP422P0Y7HoU2raasmx4r7gMGpbJQ1mEELNOAlscs2cPPMu3nvsWaxau4Wdrf4a7HB5IUii8w0tg+WfMsPh8ub9bCDErJLDFcXlwz4Pc/urtXHXyVdxy7i1FfxvEtOpvg7efMcO+FyAVNx3WTvwkLLsYlq6F8GKnaymEmKOk05k4Lteccg2dQ538etevcSs33zvne+UT2jWtcPb1ZkhEzbXud6wAf+sJs074BFiyBpaugRMulIe0CCFmhQS2KOimM28irdP8ZvdvyOgMN597c/mEdpYvCCdfYgatoWuPCfB9L5jT5689YNZrXAlLLoRF58LC86B6vrP1FkLMSRLYoiClFBvO2oBLubh/1/0MJAb40fk/wuf2OV01ZygFTaea4bx/Ms8z79gJ+56H91+A7b+GLb8w69YsssLbGppWystJhBBTJtewxYS01tz75r3c9dpdnNF4Bnd98i5qK2qdrlbxSSWg8w04uAUOvgIHtkCk0yzzhWD+amhZbY1PN+/2loe3CCFGkU5nYsqe2vcU39/8feYF5/GTNT9hZcNKp6tU3LSGvgNWgG+BQztMoKeHzXJ/NTSfZkK85XRo+oi5D9ztdbbeQghHSWCLabHzyE42vLCB7lg3G87awNUrri6/69pTkU7CkbdMeB/6u3nDWOebIyHu9kHDcmg8FRpPMafSG081HeGknYUoCxLYYtr0xfu4efPNvNj+IufNP48ffOwHtFa1Ol2t0pUN8cO7oGs3HN5tOrcNtI2s4682Ad6wDOpPgnprXLdE7g8XYo6RwBbTKqMzPLz3Ye7cficazTdO/wbrV6zH45I+jNMm1pcf5F174Og7MNQ1so5yQe0iW4ifaKbDi6FmoZxeF6IESWCLGdER6eC2V27jpfaXWFKzhO+c+R0ubL1QTpPPpHg/dL8H3e+a4eg71vR7kBwaWU+5oHqBedFJePHYcahZOr0JUYQksMWM0Vrz3MHnuHP7newf2M/ZzWdz42k3cm7zuRLcs0lrGOwwwd23H3r3548HO/LXd/tMoNe0QnWLNSzInw42SKgLMcsksMWMS6aTbHx7I/e8cQ9HY0c5reE0rl91PWta15TH88iLXTJuHrva98FIiPe3wcAhGGiHgQ7IJPO3cXnNQ2CyQR5qNk91CzVBqBEqG810ZYPcZy7ENJHAFrNmOD3MY+8+xn1v3kd7pJ35lfNZt2wdVy67knlBeYRn0cpkIHrUCu9DtiC3pvvbINKVf9o9R5nQrmw0QR5qGgn2ynkQrIdgHQTqzLS/Snq9CzEOCWwx61KZFM8eeJaH336YVzpewa3cnNdyHpeecCkXLbqIKl+V01UUH8ZwxHR8i2SHw2acV2aVZ29XG83lNQEerLdCvG4k1IP1+eUVNSOD9IgXZUACWzhq/8B+Nr2ziaf3Pc2hoUN4XV7ObzmfT7R+ggsWXEBLqMXpKorppjUMD0DkCMR6INoD0W4zxLLTtvLsOjo9/md6AvkBPnoI1I4t89eAP2SeNuerlCN7UfQksEVR0FrzxtE3ePqDp/nL/r9waOgQAEtrlvLxlo9zZtOZrG5cTUOgweGaCkdkMjDcPxLksR6ID0C8z/SOz41HDTGrfKKwB9Nz3meFt98+rjJjf9XYsuy8v9raNmhet+oNmEGu3YtpJoEtio7Wmn0D+9jctpmXDr3Ets5tJDIJABZVLWJ142pWNaxiRd0KloeXE/QGHa6xKGpaQ2JoVJj3mcBPDJpT+YmINR6E4UJl1nw6cezf66mwwtsKcV8QvNkhYI7qx10etOYDZtrjN5+XG2zz0lu/bEhgi6KXSCfY3b2bHV07eK3rNXZ07aB3uBcAhWJR9SJODp/MsvAyFlcvzg2V3kqHay7mnFTCCvLBsSGfjJkfBskYJKNmSESt+aEJlluDzny4Orl9thAP2MLcbwLfPj/ZcrffPFTH7bPmvVaZz0x7bNOj13V55LLCDJPAFiVHa83h6GHe6nmLt3reYm/PXvb07KE90p63XkOggcXVi1lYtZDmymaagk00VzbTHGymubKZkC/k0L9AiFG0htTwSHiPDvfUMKTitnE8fz4ZH395MjZqe9t8JjW9/44x4W4F+rg/BEYFvstjynLj7LTHTGeX2Ze7rXn78tx63lHbum3TBbYt8h8cEwW2PEtSFCWllAneymbWLlybK4+lYhwcPMiBgQN8MPAB+wf2c2DgAH9r/xtHYkfQ5P8ADXlDzAvOI+wPUx+op66ijnBFODeur6in1l9Lla+KKl8VQU9QHvgiZoZS4K0wA3Wz973pVH7QpxP5QypxbGV55UlzF8C46ybN5Qj7uqmEudc/kzJ1yiTNstH3/8805S7wo8FjlbtHfiy4PNZ8oTLPyOd84tvQWjBfp92UAlsp9QXgVuAU4BytdcFDYqXUJcBdgBu4R2t9x1S+V5SvgCfA8vByloeXj1mWzCQ5Ej3C4ehhOoc66Rzq5HD0MF3RLnrjvbzX9x7b4tvoG+4bE+xZLuUi5A3lArzKV5Wbr/RWUuGpIOAJEPQECXgCufnsYF/ud/vxuX34XD48Lo/8EBDOcHvAbXWeK1aZ9Eh4jwn0lFVmLU9b8wWX28a5dZPjLLd9ZiZtDamRQRcoy6QhmbDNZ8zZjVky1SPsN4ErgV+Ot4JSyg38HPg00AZsVUo9rrXePcXvFiKP1+WlJdQy6W1iqUyK/uF+euI99MR76B3uJZKIEElEGEgMEElGGEwM5oa2SBuDiUGiySixVIzkhzwi8Ll8JsDdPrwuL16XNzfvc/nwur1563hcHjzKg9vlxq3ceFwe3MqN2+XOKx89P9l6LuXKDQplxsqMXZhpt3KPlOPKWye7jVu5x2xX6HPtn5X9D8jN26ezP2py66mR6dHrizkkeyRLhdM1KWpTCmyt9R6Y9H+ec4B3tdbvW+s+BFwOSGALR3hcHuoD9dQH6j/U9qlMingqTiwVI56KE02ZII+n48SSsZHpVIxkOkkikyCZSZJIJ0ikbdMZa95aJ5FOMJQaom+4j0Q6QUqnSGVSpHWadCZNWqfz5rPLy5k9yLNhXijk7X+jxmxjrQ/kldt/PIxelvuOUX/7cnWYoGyyHxuTrX+8ywt+x2Sfebx1noZ2mKze01HnY2kb28rH5KYzbuKs5hI4JX6MFgAHbfNtwLmFVlRK3QjcCLBo0aKZr5kQH4LH5SHkCxVNh7aMzuQCfLxgt5dndIYMGbTWZlpn0IxM58q0JoNt2rZdWqfHlE34WWjSmXTepQitNbn/rOkxy2zrZBeP3sbecbbQNvZ17Nva1x+9rb0Ok61nr/No4116GW+bMZ85yXeM+fwCXzflzxz9eZPV+RjaYbJ6T7XOBf8Nx9G/erI2sJvN1wpP+k1Kqf8DmgssukVr/dh0VkZrfTdwN5he4tP52ULMVS7lwuV24UXefy3EXDZpYGutPzXF72gHFtrmW60yIYQQQhyj2Xh8zlZgmVJqiVLKB6wHHp+F7xVCCCHmjCkFtlLqCqVUG/Ax4I9KqWes8hal1JMAWusU8DXgGWAPsFFrvWtq1RZCCCHKy1R7iT8KPFqg/BDwWdv8k8CTU/kuIYQQopzJE+WFEEKIEiCBLYQQQpQACWwhhBCiBEhgCyGEECVAAlsIIYQoARLYQgghRAmQwBZCCCFKgAS2EEIIUQIksIUQQogSoAq9Cq0YKKWOAPun+WMbgKPT/JmlStoin7RHPmmPEdIW+aQ98k13eyzWWs8rtKBoA3smKKW2aa1n503jRU7aIp+0Rz5pjxHSFvmkPfLNZnvIKXEhhBCiBEhgCyGEECWg3AL7bqcrUESkLfJJe+ST9hghbZFP2iPfrLVHWV3DFkIIIUpVuR1hCyGEECWpLAJbKXWJUmqvUupdpdT3nK6P05RSHyil3lBK7VBKbXO6PrNNKXWfUqpLKfWmraxOKfVnpdQ71jjsZB1nyzhtcatSqt3aP3YopT7rZB1nk1JqoVLqOaXUbqXULqXUN63ysts/JmiLstw/lFIVSqlXlVI7rfa4zSpfopTaYuXL75RSvhmrw1w/Ja6UcgNvA58G2oCtwNVa692OVsxBSqkPgLO01mV5L6VS6kIgAvxGa/0Rq+w/gR6t9R3Wj7qw1vpfnaznbBinLW4FIlrr/3Kybk5QSs0H5mutX1NKVQHbgc8D/0iZ7R8TtMUXKcP9QymlgEqtdUQp5QU2A98Evg1s0lo/pJT6b2Cn1voXM1GHcjjCPgd4V2v9vtY6ATwEXO5wnYSDtNZ/BXpGFV8OPGBNP4D5wzTnjdMWZUtr3aG1fs2aHgT2AAsow/1jgrYoS9qIWLNea9DARcDvrfIZ3TfKIbAXAAdt822U8U5n0cCflFLblVI3Ol2ZItGkte6wpjuBJicrUwS+ppR63TplPudP/xailDoBOB3YQpnvH6PaAsp0/1BKuZVSO4Au4M/Ae0Cf1jplrTKj+VIOgS3GukBrfQZwKfDP1mlRYdHmOtHcvlY0sV8AJwKrgQ7gJ85WZ/YppULAI8C3tNYD9mXltn8UaIuy3T+01mmt9WqgFXP2dsVsfn85BHY7sNA232qVlS2tdbs17gIexex45e6wdc0ue+2uy+H6OEZrfdj6w5QBfkWZ7R/W9clHgN9qrTdZxWW5fxRqi3LfPwC01n3Ac8DHgFqllMdaNKP5Ug6BvRVYZvXk8wHrgccdrpNjlFKVVgcSlFKVwMXAmxNvVRYeB66zpq8DHnOwLo7KBpPlCspo/7A6Ft0L7NFa32lbVHb7x3htUa77h1JqnlKq1poOYDoy78EE9z9Yq83ovjHne4kDWLcd/AxwA/dprf/d4So5Rim1FHNUDeABHiy39lBK/S+wFvOWncPAD4A/ABuBRZi3xH1Raz3nO2ON0xZrMac7NfAB8FXb9ds5TSl1AfAi8AaQsYpvxly7Lav9Y4K2uJoy3D+UUqdhOpW5MQe7G7XWP7T+pj4E1AF/B76ktR6ekTqUQ2ALIYQQpa4cTokLIYQQJU8CWwghhCgBEthCCCFECZDAFkIIIUqABLYQQghRAiSwhRBCiBIggS2EEEKUAAlsIYQQogT8P76ILVv1+fTYAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 576x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def hill (S : np.ndarray, vmax : float, km : float, alpha : float) -> np.ndarray:\n",
    "    '''\n",
    "    Hill kinetic\n",
    "    \n",
    "    Parameters\n",
    "    ----------\n",
    "        S : array-like\n",
    "            Concentration of S\n",
    "        \n",
    "        vmax : float\n",
    "            Maximum velocity of the kinetic equation\n",
    "            \n",
    "        km : float\n",
    "            \n",
    "        alpha : float\n",
    "            Exponential term of the Hill equation\n",
    "    '''\n",
    "    return - (vmax * S**alpha) / (km**alpha + S**alpha)\n",
    "\n",
    "# generate a range of possible S values\n",
    "S = np.linspace(start=0, stop=30, num=1000)\n",
    "\n",
    "# classical MM solution\n",
    "dS_MM = hill(S, vmax=1, km=3, alpha=1)\n",
    "\n",
    "dS_a2 = hill(S, vmax=1, km=3, alpha=2)\n",
    "dS_a4 = hill(S, vmax=1, km=3, alpha=4)\n",
    "\n",
    "fig, ax = plt.subplots(nrows=1, ncols=1, figsize=(8, 6))\n",
    "ax.plot(S, dS_MM, label=r'$\\alpha=1$ - MM solution')\n",
    "ax.plot(S, dS_a2, label=r'$\\alpha=2$')\n",
    "ax.plot(S, dS_a4, label=r'$\\alpha=4$')\n",
    "ax.set_title('Hill Phase Space', fontsize=14)\n",
    "ax.legend()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8l43MfXzMRWF"
   },
   "source": [
    "We can easily notice how incrementing the $\\alpha$ exponent the system requires more time to achieve the critical threshold."
   ]
  }
 ],
 "metadata": {
  "colab": {
   "authorship_tag": "ABX9TyO4OCZC4wA4Dr3LUVnEbXTz",
   "collapsed_sections": [],
   "name": "lecture2_michaelis_menten.ipynb",
   "provenance": [],
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}