{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Dynamic Flux Balance Analysis (dFBA) in COBRApy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following notebook shows a simple, but slow example of implementing dFBA using COBRApy and [scipy.integrate.solve_ivp](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html). This notebook shows a static optimization approach (SOA) implementation and should not be considered production ready.\n",
    "\n",
    "The model considers only basic Michaelis-Menten limited growth on glucose."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from tqdm import tqdm\n",
    "\n",
    "from scipy.integrate import solve_ivp\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create or load a cobrapy model. Here, we use the 'textbook' e-coli core model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cobra\n",
    "from cobra.test import create_test_model\n",
    "model = create_test_model('textbook')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Set up the dynamic system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Dynamic flux balance analysis couples a dynamic system in external cellular concentrations to a pseudo-steady state metabolic model.\n",
    "\n",
    "In this notebook, we define the function `add_dynamic_bounds(model, y)` to convert the external metabolite concentrations into bounds on the boundary fluxes in the metabolic model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def add_dynamic_bounds(model, y):\n",
    "    \"\"\"Use external concentrations to bound the uptake flux of glucose.\"\"\"\n",
    "    biomass, glucose = y  # expand the boundary species\n",
    "    glucose_max_import = -10 * glucose / (5 + glucose)\n",
    "    model.reactions.EX_glc__D_e.lower_bound = glucose_max_import\n",
    "    \n",
    "\n",
    "def dynamic_system(t, y):\n",
    "    \"\"\"Calculate the time derivative of external species.\"\"\"\n",
    "\n",
    "    biomass, glucose = y  # expand the boundary species\n",
    "    \n",
    "    # Calculate the specific exchanges fluxes at the given external concentrations.\n",
    "    with model:\n",
    "        add_dynamic_bounds(model, y)\n",
    "        \n",
    "        cobra.util.add_lp_feasibility(model)\n",
    "        feasibility = cobra.util.fix_objective_as_constraint(model)\n",
    "        lex_constraints = cobra.util.add_lexicographic_constraints(\n",
    "            model, ['Biomass_Ecoli_core', 'EX_glc__D_e'], ['max', 'max'])\n",
    "    \n",
    "    # Since the calculated fluxes are specific rates, we multiply them by the\n",
    "    # biomass concentration to get the bulk exchange rates.\n",
    "    fluxes = lex_constraints.values\n",
    "    fluxes *= biomass\n",
    "    \n",
    "    # This implementation is **not** efficient, so I display the current\n",
    "    # simulation time using a progress bar.\n",
    "    if dynamic_system.pbar is not None:\n",
    "        dynamic_system.pbar.update(1)\n",
    "        dynamic_system.pbar.set_description('t = {:.3f}'.format(t))\n",
    "        \n",
    "    return fluxes\n",
    "\n",
    "dynamic_system.pbar = None\n",
    "\n",
    "\n",
    "def infeasible_event(t, y):\n",
    "    \"\"\"\n",
    "    Determine solution feasibility.\n",
    "    \n",
    "    Avoiding infeasible solutions is handled by solve_ivp's built-in event detection.\n",
    "    This function re-solves the LP to determine whether or not the solution is feasible\n",
    "    (and if not, how far it is from feasibility). When the sign of this function changes\n",
    "    from -epsilon to positive, we know the solution is no longer feasible.\n",
    "    \n",
    "    \"\"\"\n",
    "    \n",
    "    with model:\n",
    "        \n",
    "        add_dynamic_bounds(model, y)\n",
    "        \n",
    "        cobra.util.add_lp_feasibility(model)\n",
    "        feasibility = cobra.util.fix_objective_as_constraint(model)\n",
    "        \n",
    "    return feasibility - infeasible_event.epsilon\n",
    "\n",
    "infeasible_event.epsilon = 1E-6\n",
    "infeasible_event.direction = 1\n",
    "infeasible_event.terminal = True"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the dynamic FBA simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "t = 5.804: : 185it [00:16, 11.27it/s]\n"
     ]
    }
   ],
   "source": [
    "ts = np.linspace(0, 15, 100)  # Desired integration resolution and interval\n",
    "y0 = [0.1, 10]\n",
    "\n",
    "with tqdm() as pbar:\n",
    "    dynamic_system.pbar = pbar\n",
    "\n",
    "    sol = solve_ivp(\n",
    "        fun=dynamic_system,\n",
    "        events=[infeasible_event],\n",
    "        t_span=(ts.min(), ts.max()),\n",
    "        y0=y0,\n",
    "        t_eval=ts,\n",
    "        rtol=1e-6,\n",
    "        atol=1e-8,\n",
    "        method='BDF'\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Because the culture runs out of glucose, the simulation terminates early. The exact time of this 'cell death' is recorded in `sol.t_events`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  message: 'A termination event occurred.'\n",
       "     nfev: 179\n",
       "     njev: 2\n",
       "      nlu: 14\n",
       "      sol: None\n",
       "   status: 1\n",
       "  success: True\n",
       "        t: array([0.        , 0.15151515, 0.3030303 , 0.45454545, 0.60606061,\n",
       "       0.75757576, 0.90909091, 1.06060606, 1.21212121, 1.36363636,\n",
       "       1.51515152, 1.66666667, 1.81818182, 1.96969697, 2.12121212,\n",
       "       2.27272727, 2.42424242, 2.57575758, 2.72727273, 2.87878788,\n",
       "       3.03030303, 3.18181818, 3.33333333, 3.48484848, 3.63636364,\n",
       "       3.78787879, 3.93939394, 4.09090909, 4.24242424, 4.39393939,\n",
       "       4.54545455, 4.6969697 , 4.84848485, 5.        , 5.15151515,\n",
       "       5.3030303 , 5.45454545, 5.60606061, 5.75757576])\n",
       " t_events: [array([5.80191035])]\n",
       "        y: array([[ 0.1       ,  0.10897602,  0.11871674,  0.12927916,  0.14072254,\n",
       "         0.15310825,  0.16649936,  0.18095988,  0.19655403,  0.21334507,\n",
       "         0.23139394,  0.25075753,  0.27148649,  0.29362257,  0.31719545,\n",
       "         0.34221886,  0.36868605,  0.3965646 ,  0.42579062,  0.4562623 ,\n",
       "         0.48783322,  0.52030582,  0.55342574,  0.58687742,  0.62028461,\n",
       "         0.65321433,  0.685188  ,  0.71570065,  0.74425054,  0.77037369,\n",
       "         0.79368263,  0.81390289,  0.83089676,  0.84467165,  0.85535715,\n",
       "         0.8631722 ,  0.86843813,  0.8715096 ,  0.8727423 ],\n",
       "       [10.        ,  9.8947027 ,  9.78040248,  9.65642157,  9.52205334,\n",
       "         9.37656372,  9.21919615,  9.04917892,  8.86573366,  8.6680879 ,\n",
       "         8.45549026,  8.22722915,  7.98265735,  7.72122137,  7.442497  ,\n",
       "         7.14623236,  6.83239879,  6.50124888,  6.15338213,  5.78981735,\n",
       "         5.41206877,  5.02222068,  4.62299297,  4.21779303,  3.81071525,\n",
       "         3.40650104,  3.01042208,  2.6280723 ,  2.26504645,  1.92656158,\n",
       "         1.61703023,  1.33965598,  1.09616507,  0.88670502,  0.70995892,\n",
       "         0.56344028,  0.44387781,  0.34762375,  0.27100065]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot timelines of biomass and glucose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Text(0, 0.5, 'Glucose')"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAD8CAYAAADaOstiAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd0VFXXx/HvTkJCL1IE6VVQDFWqAiooCoICClhAFLH33h7ba++PHREEG1KEh6qIgCgqhF4FQk1oQXoLpOz3jzNoREgGMpM7k9mftWYx986dmZ0l8ss59xRRVYwxxhgvRXldgDHGGGNhZIwxxnMWRsYYYzxnYWSMMcZzFkbGGGM8Z2FkjDHGcxZGxhhjPGdhZIwxxnMWRsYYYzwXE8wPF5GOwDtANDBIVV8+5vWqwGCgLLATuE5Vk7P7zKioKC1UqFCQKjbGmPzp4MGDqqoh2wCRYC0HJCLRwCqgA5AMJAC9VXV5lmtGAhNUdaiIXAj0U9Xrs/vcIkWK6IEDB4JSszHG5FciclBVi3hdx4kEMyWbAYmqulZVjwDDga7HXHMWMM33fPpxXjfGGBMBghlGFYGkLMfJvnNZLQK6+Z5fCRQTkdLHfpCIDBCRuSIyNz09PSjFGmOM8Y7X/YcPAm1FZAHQFtgEZBx7kaoOVNWmqto0Jiaot7mMMcZ4IJj/sm8CKmc5ruQ79xdV3YyvZSQiRYHuqro7iDUZY4wJQcFsGSUAtUWkuojEAr2AcVkvEJEyInK0hsdwI+uMMcZEmKCFkaqmA3cC3wMrgBGqukxEnhORLr7L2gErRWQVcDrwQrDqMcYYc2IiMlhEUkRkaZZzp4nIDyKy2vdnqaB9f7jt9GpDu40x5uTlNLRbRNoA+4Fhqlrfd+5VYKeqviwijwKlVPWRYNQXOaMBfv8dfvwR2reHpk0hOtrriowxhsxMZefBI2zfd5gDh9M5lJbBoSMZHErLIPWv55kcSsvgorrlaFC5ZFDqUNWZIlLtmNNdcT1YAEOBGYCFUa78/DM8+aR7lCwJF1wAHTq4cKpVC0S8rtAYkw/tTU1j9bZ9JKbsZ8ueVFL2HSZl72FS9qWSsvcwf+4/THqmfz1U5YrF5SaMYkRkbpbjgao6MIf3nK6qW3zPt+JupwRFZHXTbd8O06bB1Knwww+wYYM7X6WKC6ajj9NOC1zBxpiIsP9wOqu27WP1tn2s2rbf93w/W/em/uO604rEUq5YHOWKF3R/+h5lixWkWMEYCsVGU6hANAULRP/1vFCBaOJiooiKOvVfmv1ZgcHXMpqQpZtut6qWzPL6LlUNyn2jyAqjrFRhzRoXSlOnupDavRuioqBVK7jsMveIj7dWkzHmX/YcSiNh3U5+X7uD39ftYNnmvRz957RggShqlStKnXLFqH16MeqcXpTa5YpRvkRBYmO8md55imG0EminqltEpAIwQ1XPDEp9ERtGx8rIgDlzYPJkmDgR5s935ytW/DuYLroIihUL/HcbY0LeicInNiaKJlVK0az6adSvWII6pxelUqnCROeiFRMMpxhGrwE7sgxgOE1VHw5KfRZGJ7BlC3z3HUyaBFOmwN69EBvr7jFdcQV06QKnB6371BgTAvYcTOP75VsZv2gzv67ZQUam/hU+LWqUpkWN02hQuSQFC4T+gCg/RtN9jRusUAbYBjwNjAVGAFWADcDVqrozKPVZGPkhLQ1mzYLx42HMGFi3znXdtWoFV17pwqlmzbytyRgTFPsPpzN1+TbGL9rMzNXbSctQqpxWmE7xFWhXp2zYhM+xQn3Vbgujk6UKS5bA2LEumBYudOfr13fB1LMnnH22d/UZY05aWkYmPyzfxriFm5m+MoXD6ZmcUaIgneIrcHmDMzinYgkkzO8dWxgFmOdhdKz1610wjR3rho9nZrow6tnTPerU8bpCY8wJ7D54hK/mbGTYrxvYujeVMkXj6Bxfgc7xFWhcpVSuRq+FGgujAAu5MMpq61YYPRq++cYFE0DDhtCrF1x9NVSv7m19xhgAElP2M2TWOkbPTyY1LZPzapWhX+tqtDuzXMgNPAgUC6MAC+kwyio5GUaOdME0e7Y717w5XH+9C6fS/9q2yRgTRKrKL4l/8ukv65ixcjuxMVFc2bAi/c6rRt3yxb0uL+gsjAIsbMIoq/XrYcQI+PJLWLwYChSAzp2hb1+49FI3Ss8YExSqypTl23hzyipWbttHmaJx9GlZlWuaV6FM0Tivy8szFkYBFpZhlNWiRTBsmAumbdugTBno3dsFU+PGNsHWmACau34nL03+g3kbdlGzbBFua1eLyxtUIC4m/EbD5ZaFUYCFfRgdlZ7u5i8NHQr/+x8cPuwGPvTvD3362JJExuRCYso+XvluJT8s30a5YnHc16EOVzWpREy015tbe8fCKMDyTRhltXu368YbPNjdX4qLg6uuggED4LzzrLVkjJ+27U3l7amr+CYhicKxMdzatgY3nledwrGRsyb0iVgYBVi+DKOsFi2CTz6Bzz93qz7UretCqU8fG/RgzAkcOpLBBzMS+eTntWRkKtc2r8pdF9aidATdE8pJRIeRiHQE3gGigUGq+vIxr1fB7ZFR0nfNo6o6KbvPzPdhdNSBA661NHCg24spLg569IC77nKj8owxAPya+CePfruEjTsPcnmDM3jw4jpULR2y/+Z6JmLDSESigVVAByAZSAB6q+ryLNcMBBao6ocichYwSVWrZfe5ERNGWS1e7FpLw4a51lKzZnDPPS6cbCSeiVB7DqXx0qQVDE9IolrpwrzULZ6WNa334ERCPYyCeTevGZCoqmtV9QgwHLdrYFYKHB3gXwLYHMR6wld8PLz7rpu79N577h7TtddC1arw3HNuVJ4xEeS7pVvp8OZPjJibxC1ta/DdvW0siMJcMFtGPYCOqtrfd3w90FxV78xyTQVgClAKKAK0V9V52X1uRLaMjpWZ6Ubi/fe/bsuL2Fg3kfaee9zwcGPyqZR9qTwzbhmTlmylXoXivNo9nnMqlfC6rLAQyS0jf/QGPlPVSsBlwOci8q+aRGSAiMwVkbnp6el5XmTIiYqCjh3d9hZ//OEGOIweDU2auC0ufvgBwmxgijHZUVVGzk2iw5szmboihYcuOZNxd7a2IMpHgtkyagk8o6qX+I4fA1DVl7JcswzXekryHa8FWqhqyok+11pGJ7Bnj7uv9NZbsHkzNGoEjzwC3btDjA1rNeFrb2oaj45ezKQlWzm3Wile7h5PzbJFvS4r7ERyyygBqC0i1UUkFugFjDvmmo3ARQAiUg8oCGwPYk35V4kS8OCDsHYtDBoEBw+6rrszz4QPP4RDh7yu0JiTtihpN53++zPfL9vGo5fW5ZsBLS2I8qlgD+2+DHgbN2x7sKq+ICLPAXNVdZxvBN0nQFHcYIaHVXVKdp9pLSM/ZWa6lR1eecVNpC1Xzt1TuvNOKJ7/F4U04U1VGTxrPS9PXkHZonG8e00jmlS1VUlyI9RbRjbpNb9ThZkzXShNngylSsH998Pdd1somZC0++ARHhy5mKkrttG+3um8flU8JQvbFIbcsjAKMAujXJg3D5591m2fXqoUPPCAm0RroWRCxLwNu7j76wWk7Evl0UvrcWPramG/w2qosDAKMAujAJg3D555BiZMcAuyHg2lYsW8rsxEKFVl4My1vPb9SiqULMh7vRvToHJJr8vKVyyMAszCKIDmznWhNHGiC6UHH3Tdd0VC9u+ryYdS0zJ4cOQiJizewqX1y/Ny93hKFCrgdVn5joVRgFkYBUFCguu+mzgRypd3AXXjjW4TQGOCaNveVAYMm8viTXt4+JK63Nq2hnXLBUmoh5HXk15NKDj3XNdl98svULMm3Hor1K8Po0bZ5FkTNEs37aHre7NYnbKfj69rwm3taloQRTALI/O31q3h55/dkPCYGLenUosWMGOG15WZfGbSki30+OhXoqOEUbe24uKzy3tdkvGYhZH5JxHo0sXtq/Tpp7BpE1xwAVx2mVs93JhcUFX+++Nqbv9yPmdVKM7YO1pz1hk2mtNYGJkTiYlx941Wr3ZzlH77DRo2dOvgbbdFMszJS03L4O7hC3nzh1Vc2agiX93cgrLFbPM749gABuOfXbvcdhXvvedG2z39NNxxh+2nZPzy5/7D3DR0LouSdvPQJWdyu90fynOhPoDBwsicnBUr4L774Pvv3bp3b7/tVhA35gSSdh7k+k9ns3VvKm/3bETH+nZ/yAuhHkbWTWdOTr16blmhCRPc+neXXgqdO8OqVV5XZkLQii176fbhr+w6mMaX/ZtbEJkTsjAyJ08EOnWCpUvh9dfd2nf167tJs/v2eV2dCRGz1+7g6o9/I1qEkbe2tIVOTbYsjMypi411SwmtXg19+sCbb8JZZ8GYMTY/KcJNWbaV6wfPoWyxOEbf3oo6p9tSUyZ7FkYm904/3e2h9Ouvblmhbt2ga1fYsMHryowHRiQkcesX86hXoTijbm1FxZKFvC7J+EFE7hORZSKyVES+FpGCefn9FkYmcFq0cIuwvv46/PijayW99hqkpXldmckDqsoHMxJ5ePRizqtdlq/6N+e0IjbaMhyISEXgbqCpqtbH7UHXKy9rsDAygRUT47ruVqyA9u3h4YehSRM3T8nkW5mZyv9NXMGr362kS4MzGNSnKUXibLv7MBMDFBKRGKAwsDkvv9zCyARHlSpuWaExY9wcpdat3Zp3e/Z4XZkJsMxM5fExS/j0l3Xc0Koab/dsSGyM/dMSTlR1E/A6sBHYAuzJadftQAvq3xgR6SgiK0UkUUQePc7rb4nIQt9jlYjsDmY9xgNXXPH33KRPPnGj7r77zuuqTIBkZCoPjlrE8IQk7rqwFk9ffhZRUTaZNUTFiMjcLI8BR18QkVJAV6A6cAZQRESuy8vigjbpVUSigVVAByAZSAB6q+ryE1x/F9BIVW/M7nNt0msYmzMH+vWD5cvdn2++CSVtA7VwlZ6RyX0jFjF+0Wbu71CHuy+q7XVJJhvZTXoVkauAjqp6k++4D9BCVW/Pq/qC2TJqBiSq6lpVPQIMxyXvifQGvg5iPcZrzZrB/Pnw+OMwbBicfbbbQ8mEnSPpmdz19QLGL9rMo5fWtSAKfxuBFiJSWNw6TRcBK/KygGCGUUUgKctxsu/cv4hIVVzzcNoJXh9wtGmZnp4e8EJNHoqLgxdegNmz3TDwzp3hhhvcfSUTFg6nZ3D7l/OYvHQrT3U+i1vb1vS6JJNLqjobGAXMB5bgsmFgXtYQKncZewGjVDXjeC+q6kBVbaqqTWNibIROvtCkidv2/Mkn4YsvXCtp/HivqzI5SE3LYMCweUxdkcLzXc/mpvOqe12SCRBVfVpV66pqfVW9XlUP5+X3BzOMNgGVsxxX8p07nl5YF13kiYuD559395LKlnX7KPXvD/v3e12ZOY5DRzLoP3QuM1dv5+Vu53B9y2pel2TykWCGUQJQW0Sqi0gsLnDGHXuRiNQFSgE2ESVSNW4MCQnw2GMweDA0auS68UzIOHA4nRuGzOHXNX/yeo8G9GpWxeuSTD4TtDBS1XTgTuB73I2wEaq6TESeE5EuWS7tBQzXcNvLwgRWbCy8+KLb4vzIETcv6fnnwe4Reu7QkQxuGppAwvqdvNWzId2bVPK6JJMP2X5GJvTs3u027vvqK2jVyt1Tqm73JryQmua65n5d8ydv9WxI14bHHYNkwoDtZ2TMySpZEr780j2WLoUGDdxQ8DD7xSncpaZlcMvn85i15k9e69HAgsgElYWRCV3XXAOLF0PDhtC3L/Tq5VpNJuiOpGdy+5fz+WmVG6xgXXMm2CyMTGirWhWmT3f3k7791g0JnzfP66rytbSMTO78aj7T/kjh/66oT89zbbCCCT4LIxP6oqPdSLuZM912FK1awXvvWbddEKRnZHLP8AVMWb6NZ7uczXUtqnpdkokQFkYmfLRsCQsWQIcOcNddcPXVtgp4AGVkKvePWMSkJVt5slM9+raq5nVJJoJYGJnwUro0jBsHr77qtqdo3Nitd2dyJSNTeWjkIsb51prrf34Nr0syEcbCyISfqCh46CHXbXfkiGsxffCBddudIlXlybFL+HbBJh7oUMfWmjOesDAy4atVK1i40O0oe8cd0LMn7NvndVVhRdXt0Pr1nCTuuKAmd9nq28YjFkYmvJUu7RZYffllN9queXNYudLrqsLGW1NX/7VD64MXn+l1OSaCWRiZ8BcVBY88Aj/8ANu3u32TbAXwHH380xr+++Nqrm5aif90Pgu3jY0x3rAwMvnHBRe4OUi1a7sVwJ99FjIzva4qJH3++wZemvwHneMr8FK3eNsq3HjOwsjkL1WqwM8/Q58+8MwzcMUVNvz7GKPnJfPU2KW0r1eOt3o2JNqCyIQACyOT/xQqBJ99Bu++C5Mnu2675cu9riokTF6yhYdGLaJ1rdK8d01jCkTbPwEmNNjfRJM/icCdd8K0aW49u+bN3QCHCDZ9ZQp3D19AoyqlGHh9UwoWiPa6JGP+YmFk8rfzz3f3kc46C7p3h+eei8j5SLPX7uDWz+dR5/RiDL7hXIrExXhdkjH/YGFk8r9KldwE2b594emn4dpr4dAhr6vKM0uS93DT0LlUKlWIYTc2o0ShAl6XZMy/BDWMRKSjiKwUkUQRefQE11wtIstFZJmIfBXMekwEi4uDIUPgpZfg66/dyLutW72uKuhWb9tHn8GzKVm4AF/2b0HponFel2TMcQVtp1cRiQZWAR2AZCAB6K2qy7NcUxsYAVyoqrtEpJyqpmT3ubbTq8m1MWPguuv+njDboIHXFQVF0s6D9PjoVzIVRt7SkmplQnaTT5MHInmn12ZAoqquVdUjwHCg6zHX3Ay8r6q7AHIKImMC4sor4Zdf3Byk1q3z5QTZlL2pXPfpbFLTMvn8pmYWRCbkBTOMKgJJWY6TfeeyqgPUEZFZIvK7iHQMYj3G/K1RI5gzB+rVg65d4fXX883Aht0Hj3D9p3PYvu8wQ/qdS93yxb0uyZgceT2AIQaoDbQDegOfiEjJYy8SkQEiMldE5qanp+dxiSbfOuMM+Okn6NHDrQLev79bBTyM7T+cTt8hCaz78wCf9GlK4yqlvC7JGL8EM4w2AZWzHFfyncsqGRinqmmqug53j+lfywar6kBVbaqqTWNibEiqCaDChWH4cHjqKRg8GDp1gr17va7qlKSmZTBg2FyWbtrDu9c0onWtMl6XZIzfghlGCUBtEakuIrFAL2DcMdeMxbWKEJEyuG67tUGsyZh/i4py84+GDIEZM9zcpE3H/t4U2tIyMrnr6wX8umYHr/WI55Kzy3tdkjEnJWhhpKrpwJ3A98AKYISqLhOR50Ski++y74EdIrIcmA48pKo7glWTMdm64QaYMAHWrnUb9i1b5nVFfsnMVB4etZgflm/j2S5n061xJa9LMuakBW1od7DY0G4TdAsWwGWXuYmxY8dCu3ZeV3RCqsrT45Yx7LcNPNChjm2OZ04okod2GxOeGjWC3393AxwuucTdUwpRb/6wimG/beDm86tz54W1vC7HmFNmYWTM8VStCrNmQYsW0Lt3SA79/mTmWt6dlkjPppV5/LJ6tjmeCWsWRsacSKlS8P33cPXVbuj3PfdARobXVQEwfM5GXpi0gk7nVODFbudYEJmwZ+OkjclOwYJuLbvKleGNN2DbNvj8c4iN9aykiYu38NiYJbStU9Y2xzMB5ZvnOQioDyhwo6r+lhffbWFkTE6iolw3XfnyroW0ezeMHg1Fi+Z5KTNWpnDvNwtoUqUUH13XhNgY69wwAfUO8J2q9vBNySmcV1/s199kEWqKEOd73k6Eu0X410oJxuRrDz7oJsZOnQrt28OOvJ2FkLB+J7d+MY/a5Yrx6Q3nUijWNsczgSMiJYA2wKcAqnpEVXfn1ff7+2vVaCBDhFrAQNzKCrbdg4k8/fq5VtHChdCmTZ5Njl26aQ83fpbAGSUKMewm25PIBEV1YDswREQWiMggETn5oeAiVRFp73teCJFi/rzN3zDKVCUduBJ4V5WHgAonXaQx+cEVV8DkyZCU5Fb9Xr06qF+XmLKPPoPnULxgAT7v35wytieROTUxR9f49D0GHPs60Bj4UFUbAQeA4+5Dd0IiNwOjgI99ZyrhVtrJkb9hlCZCb6AvMMF3zn41M5Hrggtg+nQ4cMAF0oIFQfmapJ0HuW7QHKJE+KJ/cyqWLBSU7zERIf3oGp++x8BjXk8GklV1tu94FC6cTsYdQGvALfCouhoo588b/Q2jfkBL4AVV1olQHfj8JIs0Jn9p0sTti1SoELRt61YAD6Bte1O5dtBsDqVl8EX/ZlS3PYlMEKnqViBJRM70nboIWJ7NW47nMG7/OkckBjcqL0cnvRyQCKWAyqosPqk3BogtB2RCTnIyXHyxW9Nu1Cjo3DnXH7nzwBF6fvwbm3cf4subW9Cwso0XMrnjz3JAItIQN7Q7Frdodb+jm5/6+SWvAruBPsBdwO3AclSfyPGt/oSRCDOALrg+xXlACjBLlfv9LjJALIxMSNqxAzp2dAMbvv7a7ZF0ivampnHtJ7NZtW0fn/VrRsuapQNYqIlUebI2nUgUcBNwMSC4xbAH4UfQ+NtNV0KVvUA3YJgqzYH2p1iuMflP6dJuyHfz5tCzp5sYewoOHcngps8SWLFlLx9d18SCyIQX1UxUP0H1KmAAMNufIAL/wyhGhArA1fw9gMEYk1WJEm75oHbtoG9f+PjjHN+S1eH0DAZ8Ppd5G3bxTq9GXFDXr/u+xoQOkRmIFEfkNFwv2ieIvOXPW/0No+dwza1EVRJEqAEEdzyrMeGoSBG3J9Kll8Ktt8Jbfv1/SHpGJvd8vZCfV//Jy93i6RRvMydMWCqB6l+9aKg2xw2EyJFfYaTKSFXiVbndd7xWle6nXK4x+VmhQjBmDHTvDvffDy+8kO3lGZnKQ6MW892yrfyn81lcfW7lPCrUmICLQeSUetH8WptOhIK4m1JnAwWPnlflxpP5MmMiRmys2wepXz948kk3H+mFF+CY1bUzM5UnxixhzIJNPHTJmdx4XnWPCjYmII72os1CNQERv3vR/O2m+xwoD1wC/ISbVbsvpzeJSEcRWSkiiSLyr5m8InKDiGwXkYW+R38/6zEm9MXEwNChcPPN8NJLcO+9/9gTSVV5dvwyhickcfeFtbjjAtscz4Q51ZGoxqN6m+94Lap+9aL5u2p3LVWuEqGrKkNF+Ar4Obs3iEg08D7QATezN0FExqnqsZOovlHVO/2sw5jwEhXlBjIULgzvvANpafDee6gIL0/+g6G+XVrv61DH60qNyT2RSsC7uFUYwOXEPagm5/RWf8MozffnbhHqA1vJeYmHZkCiqq51NcpwoCsnP6PXmPAm4gYyFCwIr7wCGRm81e0+Pp65lj4tq9ourSY/GYJbRPsq3/F1vnMdcnqjv2E00LfywlPAOKAo8J8c3lMRSMpynAw0P8513UWkDbAKuE9Vk45zjTHhTcR11UVHw4svUuH39fR64hWeufxsCyKTn5RFdUiW488QudefN/o7mm6QKrtU+UmVGqqUU+WjUyr1n8YD1VQ1HvgBGHq8i0RkwNGVZtPT0wPwtcZ4QIRPO97EO6160XvxFF6c/A5Rmul1VcYE0g5ErkMk2ve4DvBr4y9/R9OVxK01VC3re1S5O5u3bcLte3RUJd+5v6hq1iIHAa8e74N8q8sOBLcckD81GxNqvpy9gecnruDS2x4ko/2ZRD/3LGRmug37om2jPJMv3Ii7Z/QWboHUX3ELbefI3266ScDvwBLA31/lEoDaIlIdF0K9gGuyXiAiFVR1i++wC7DCz882JqyMmJvEE2OWclHdcrzTqxHRMU0gJhr+8x8XSJ99ZoFkwp/qBty/5SfN3zAqeLKLoqpquojciRtzHg0MVtVlIvIcMFdVxwF3i0gXIB3YCdxwMt9hTDgYOTeJR0Yv5vzaZXj/2sbExvh6x596ygXQE0+4QBo61A0HNyZciQzFjZ7b7TsuBbyBao5zUv1dtfs+YD9uRu3ho+dV2XmKJZ8yW7XbhJPR85J5cNQizqtVhk/6NKVggeO0fl55BR59FHr1cgusWiCZIMijVbsX4HaJzf7ccfj7t/4I8BrwBH9vlKRAjZMo05iIMmaBC6JWNUufOIgAHnnEtZAeesgdWyCZ8BWFSCmO7oHkFkz16y+zv3/jH8BNfP3z1OozJrL8b+EmHhixiBbVSzOoz7knDqKjHnzQrc7w8MNuouywYXYPyYSjN4DfEBmJ28+oB5D94ow+/oZRInDw1GozJrKMW7SZ+75ZSLPqp/HpDU0pFOtnqDz0EGRkwGOPuSAaMsQCyYQX1WGIzAUu9J3pxr9X3Tkuf8PoALBQhOn8855RdkO7jYk44xdt5t7hC2ha7TQG33AuhWNPsrvt0UddID35pAuiTz91LSVjwoFIFdz4gnH/OKe6Mae3+vt/yljfwxhzAhMXb+HebxbStOppDDmVIDrqiScgPR2eecYF0sCBFkgmXEzk73EFhYDqwErcjg/Z8uv/Ft/iqLHA0dUcV6r+tV6dMRFvwuLN3DN8IY0ql2Rwv3MpEpfLAQhPP+1aSM8/7wLpww8tkEzoUz3nH8cijcHtg5cTf1dgaIdbqmc97qZUZRH6qjLzZOo0Jj8asyCZB0YsoknVUgzp14yiuQ2io5591gXSiy+6QHr//X/th2RMSFOdj8jx1iT9F3//r3kDuFiVlQAi1AG+BpqcWoXG5A8jEpJ45NvFbtRc36a5bxFlJQL/93+uy+7VV13L6N13LZBM6BLJujhCFNAY2OzPW/39P6fA0SACUGWVCAX8r9CY/OeL3zfw5NilnF+7DAOvP4lRcydDBF5+2bWQ3njD7SD7xhsWSCZUFcvyPB13D2m0P2/0N4zmijAI+MJ3fC0w1+/yjMlnhsxax7Pjl3Nh3XJ8cG3jnOcR5YYIvPaaayG99RbExbmuOwskE2pUnz3Vt/obRrcBd8BfQ7l/Bj441S81JpwNnLmGFyf9wSVnn867vbOsNRdMRzfoO3zYtZTi4txoO2NCgch4/h5F92+qOS6e6u9ousPAm76HMRHrvWmreX3KKjrHV+Ctng0pEJ2HI9xE3CCGI0dS2jdEAAAah0lEQVTc4IbYWHj88bz7fmNO7PXjnDsaTn414bMNIxFGqHK1CEs4TuqpEu/PlxgT7lSVt6au5r8/rqZbo4q82iOemLwMoqOioty8oyNH3HykuDh44IG8r8OYfyoJVEL1fQBE5gBlcbnxiD8fkFPL6B7fn51PsUBjwp6q8sLEFQz6ZR1XN63ES93iiY7y8H7N0aWCDh92a9rFxsJdd3lXjzHwMG7PuqNigaZAEWAIMDKnD8g2jFTZ4vtzw9FzIpQBdqhm0z9oTD6RnpHJ42OWMGJuMn1bVuXpy88myssgOiomBr78EtLS4O67XQtpwACvqzKRKxbVpCzHv+B28t6BiF/bVmTbzyBCCxFmiPCtCI1EWAosBbaJ0PHU6zYm9B1Oz+DOrxYwYm4yd19Um2e6hEgQHVWgAAwfDpddBrfc4naLNcYbpf5xpHpnlqOy/nxATp3e7wEv4ia4TgP6q1IeaAO85HeZxoSZg0fS6T90Lt8t28pTnc/i/g51kFAcSh0XB6NHQ4cOcOON8NVXXldkItNsRG7+11mRW4A5/nxAtju9irBQlYa+5ytUqZfltQWqZLt7n4h0BN7BbTs+SFVfPsF13YFRwLmqmu38Jdvp1QTbnoNp9PtsDguTdvNK93iualrZ65JydvAgdOoEP/8M33wD3bt7XZEJMf7s9Coi0bg5pJtU1f+xAiLlcItpHwbm+842AeKAK1DdltNH5NQyyszy/NAxr2V7z8j3Q70PXAqcBfQWkbOOc10x3ECJ2TkVa0ywpexLpefA31i6aS8fXNs4PIIIoHBhGD8emjd325ePH+91RSY83QOsOOl3qaag2gp4HreG6XrgOVRb+hNEkHMYNRBhrwj7gHjf86PH5+Tw3mZAoqquVdUjwHCg63Guex54BUj1p2BjgiVp50Gu+ug3Nu48yOAbzqVj/Qpel3RyihaFSZOgYUPo0QOmTPG6IhNGRKQS0AkYdMofojoN1Xd9j2kn89Zsw0iVaFWKq1JMlRjf86PHOa1NVxHIOroi2XfuL+KWF6+sqhOz+yARGSAic0Vkbnp6eg5fa8zJW7VtHz0++pXdB9P4on9zzqtdxuuSTk2JEvD991CvHnTtCjNmeF2RCR9v44ZoZ+Z0YTB4tkGKiEThVnTIccaeqg5U1aaq2jQmJoCrIhsD/L52Bz0+/BVV+OaWFjSuUirnN4Wy006DH36AGjWgc2eYNcvrikxoiDn6S73v8ddcABHpDKSo6jyvigtmGG0Csna4V/KdO6oYUB+YISLrgRbAOBFpGsSajPmHCYs30+fTOZQrXpBvb29F3fLFvS4pMMqWhR9/hIoV4dJLYY5fA5pM/pZ+9Jd632NgltdaA118/xYPBy4UkS+O+ylBku1oulx9sEgMsAq4CBdCCcA1qrrsBNfPAB600XQmrwz+ZR3PT1xOkyqlGNS3KSULx3pdUuAlJ0ObNrBrF0yf7u4nmYjkz2g633XtcP8W5+nKO0FrGalqOnAn8D1udMYIVV0mIs+JSI4ruBoTLJmZyouTVvDchOVcclZ5vujfPH8GEUClSjBtGhQrBu3bw9KlXldkzHEFrWUULNYyMrlxOD2Dh0YuZtyizfTxLe/j6TpzeSUx0bWQMjLgp5+gbl2vKzJ5zN+WkVc8G8BgTF7bm5pGvyEJjFu0mUc61uXZLhESRAC1arkWkghceCGsXu11Rcb8g4WRiQhb9hzi6o9+Y866nbx5dQNua1czNJf3Caa6dWHqVLe46oUXwrp1XldkzF8sjEy+tyR5D1e8P4uknQcZ0u9cujWu5HVJ3qlf3wXSgQNwwQWwcaPXFRkDWBiZfO67pVu46uNfiYmKYvTtrTi/tl8LCOdvDRq4eUi7d7sW0qZNOb/HmCCzMDL5kqrywYxEbv1iPvUqFGfsHa3zzxyiQGjSxK3UkJLiAmnrVq8rMhHOwsjkO0fSM3lo1GJe/W4llzc4g69vbkHZYnFelxV6mjd3a9lt2gQXXQTbt3tdkYlgFkYmX9l14AjXfTqbUfOSueei2vy3V0MKFoj2uqzQdd55MGGCG8zQvj3s2OF1RSZC2Twjk2+s2b6fGz9LYMueVF7rEU/XhhVzfpNxpk5169jVq+eely7tdUUmwGyekTF5YOaq7Vz5/iz2p6bz9c3NLYhOVvv2MG4crFjhdo3dudPrikyEsTAyYU1V+einNdwwZA4VShRi7B2taVL1NK/LCk8XXwz/+x8sX+7CyQLJ5CHrpjNh6+CRdB4atZiJi7fQKb4Cr3aPp0icbTGSa9995/ZCql/fDQE/zcI9Pwj1bjoLIxOWNuw4wC2fz2PVtn083LEut7SpEXkrKgTT5MlwxRV/T5ItFeZ7PBkLo0CzMDI/rdrO3V8vAODd3o1oU8cmsgbFpElw5ZVwzjkukEqW9LoikwuhHkZ2z8iEDVXlwxlr6DdkDhVKFGT8nedZEAXTZZfBt9/CkiVuUMPu3V5XZPIxaxmZsHDgcDoPj1rMxCXu/tBrPeIpHGv3h/LEhAnQrZvbmG/KFGshhalQbxlZGJmQ98fWvdz+5XzW/3mARzrWZYDdH8p748dD9+6uy27KFJuHFIZCPYyC2k0nIh1FZKWIJIrIo8d5/VYRWSIiC0XkFxE5K5j1mPCiqoxISOKK92exLzWdL/o355a2Ebj1Qyi4/HIYOxaWLXOrfaekeF2RyWeC1jISkWhgFdABSAYSgN6qujzLNcVVda/veRfgdlXtmN3nWssoMhw8ks6TY5fy7fxNtK5Vmrd7NrL15ULB1KnQpQtUqwY//ggVKnhdkfFTJLeMmgGJqrpWVY8Aw4GuWS84GkQ+RYDw6jM0QbFq2z66vDeLMQs2cW/72gy7sbkFUaho397NQ0pKctuYJyV5XZHJJ4IZRhWBrH9Tk33n/kFE7hCRNcCrwN1BrMeEgVHzkun63ix2HzzCFzc15972dSJna/Bw0aaNu2+UkuKe246xJgA8H9qtqu+rak3gEeDJ410jIgNEZK6IzE1PT8/bAk2eOHgknYdGLuLBkYtoULkEk+4+n9a1ynhdljmRli1dN92ePdC2Laxe7XVFJswF855RS+AZVb3Ed/wYgKq+dILro4Bdqloiu8+1e0b5z5LkPdwzfAHrdhzgrgtqcY+1hsLHwoVuDlKBAi6c6tXzuiJzApF8zygBqC0i1UUkFugFjMt6gYjUznLYCbBfryJIRqby/vRErvxgFofSMviyf3Puv/hMC6Jw0rAhzJgBmZmuhbRokdcVmTAVtFmDqpouIncC3wPRwGBVXSYizwFzVXUccKeItAfSgF1A32DVY0JL8q6D3D9iEXPW7aRTfAVevOIcShQu4HVZ5lScfTbMnOl2i23b1s1JOv98r6syYcYmvZo897+Fm3hy7FJU4dkuZ9OtcUWbO5QfbNzotqHYsAFGjnSb9ZmQEerddBZGJs/sOZTGU2OXMm7RZppWLcVbPRtS+bTCXpdlAmn7drem3YIFMHgw9OnjdUXGJ9TDyBb3Mnni18Q/eWjUYrbuTeWBDnW4rV1NYqI9H8xpAq1sWZg2zW0/0bev26Dv3nu9rsqEAQsjE1T7UtN4afIffDV7I9XLFGHUrS1pVMX2xsnXihVz209ccw3cdx/8+Sc8/zxYV6zJhoWRCZqZq7bz2LdL2LznEDefX537O5xJodhor8syeSEuDkaMgNtugxdecIH0/vsQbf/9zfFZGJmA25uaxgsTVvDN3CRqlC3CqFtb0aSqtYYiTnQ0fPwxlCkDL70EO3bAF1+4oDLmGBZGJqCmr0zh8W+XsG1vKre0rcF97etQsID9NhyxRODFF10gPfAAbNsGY8bYFhTmX2w0nQmI3QeP8H8TVzBqXjK1yxXl1R7xdm/I/NPw4XDDDVClCkycCLVr5/gWEzihPprOwsjkiqoyZsEmXpi4gt2H0rilTQ3uvqi2tYbM8c2aBV27gqrbH8kmx+aZnMJIRCoDw4DTcTsoDFTVd/KsPgsjc6rWbN/Pk2OW8tvaHTSoXJIXr6zP2Wdku7SgMZCYCJ06wfr1MGSIG3Vngs6PMKoAVFDV+SJSDJgHXJF1D7pgsntG5qSlpmXwwfREPvppLXEFovi/K+rTu1kVW1PO+KdWLfjtN7jySrj2WlizBp580oZ+e0xVtwBbfM/3icgK3LY/FkYm9Py8ejtPjV3K+h0H6drwDJ7oVI9yxQp6XZYJN6ed5vZEuvlm+M9/XGvpk08gNtbrygwgItWARsDsvPpOCyPjl217U3lh4grGLdpM9TJF+OKm5pxX2/YbMrkQFwdDh7qW0tNPuzXtRo+2kXbBEyMic7McD1TVgcdeJCJFgdHAvcfsxh1Uds/IZCs1LYNPf1nH+9MTSc9QbmtXk9va1bQBCiawvvwSbrwRKlSAb7+Fxo29rijf8Wc0nYgUACYA36vqm3lTme+7LYzM8agq3y3dyguTVpC86xCXnH06j19Wj6qlQ3ZkqAl3c+ZA9+5usdUPP4R+/byuKF/xYwCDAEOBnaqa5wsKWhiZf1m+eS/PTVjG72t3cubpxfjP5WfZFuAmb2zfDr16ucVWb7kF3nnHVmwIED/C6DzgZ2AJkOk7/biqTsqT+iyMzFE79h/mjR9WMXzORkoUKsD9F59J73Mr2+raJm+lp8MTT8Crr0Lz5jBqFFSq5HVVYc8mvQaYhVHgpaZl8Nmv63l/eiKHjmRwfcuq3HtRHdt51Xhr9Gi3YkOhQvDNN3DBBV5XFNYiOoxEpCPwDm7b8UGq+vIxr98P9AfSge3Ajaq6IbvPtDAKnPSMTEbNS+btqavZujeVC84syxOd6lGrXDGvSzPGWbHCzUdKTISXX3br29l8pFMSsWEkItHAKqADkAwkAL2zzuYVkQuA2ap6UERuA9qpas/sPtfCKPdUle+XbeO17/9gzfYDNKpSkkc71qV5DRtSa0LQ3r1uMMO337qVGz79FE4/3euqwk6oh1EwbwY0AxJVda2qHgGGA12zXqCq01X1oO/wd8A6hoPs97U76Pbhr9z6xTwAPrquCd/e1sqCyISu4sXdfaN33oGpU+Gcc2D8eK+rMgEWzDCqCCRlOU72nTuRm4DJx3tBRAaIyFwRmZuenh7AEiPH0k176DdkDr0G/s6W3am80v0cvr+3DR3rl0es28OEOhG4+26YNw/OOAO6dIFbbwXrJck3QmIFBhG5DmgKtD3e675ZwgPBddPlYWlhb0nyHt75cRVTV6RQvGAMj15alxtaVbNJqyY8nX02zJ4NTz0Fr78O06e7CbNNm3pdmcmlYIbRJqByluNKvnP/ICLtgSeAtqp6OIj1RJTFybt5Z+pqfvwjhRKFCvBAhzr0bV2N4gVthJwJc3Fxbtj3pZdCnz7QsiU88ww8+qhtax7GgjmAIQY3gOEiXAglANeo6rIs1zQCRgEdVXW1P59rAxiytyhpN+/8uJppvhC6+fzq9G1VjWIWQiY/2rULbr/dbdzXujUMHgx16nhdVUgK9QEMwR7afRnwNm5o92BVfUFEngPmquo4EZkKnINv2XJgo6p2ye4zLYyOb96GXbw/PZFpf6RQsnABbj6/Bn1aVrUQMpHhq69cKB06BI895lpJBW01+awiOoyCwcLob5mZyrQ/Uvh45hoS1u/6K4T6tqpG0biQuB1oTN7ZutXNQ/rqK7el+QcfQPv2XlcVMiyMAszCCA6nZ/C/hZsZOHMtiSn7qViyEDefX52rz61M4VgLIRPhpk51raTVq90usm+8AeXLe12V5yyMAiySw2hfahpfzd7I4Fnr2Lb3MPUqFOfWtjXodE4FWz/OmKxSU92KDS+95JYTeuklGDAgogc4WBgFWCSG0cYdBxn223q+SUhi3+F0WtcqzS1tanJ+7TI2R8iY7Kxa5VpJP/4IzZrB22+70XcRyMIowCIljFSVXxL/ZOiv6/nxjxSiRehYvzwD2tQgvlJJr8szJnyowtdfu/tJW7dC167wwgtuzlIEsTAKsPweRgcOp/Pt/GSG/raBxJT9lC4SyzXNq3Bt86qUL2Gjg4w5ZQcOuCWFXnkF9u1zc5SefRaqVvW6sjxhYRRg+TWM1mzfz5e/b2TkvCT2paYTX6kEN7SqRqf4CsTFRG4/tzEBt2OHu5/07ruu1XTbbfD441CunNeVBZWFUYDlpzBKTcvgu6Vb+WrORuas20lMlHDZORW4oXU1GlUuafeDjAmm5GTXMho8GAoXdt1499wDpUp5XVlQWBgFWH4Io5Vb9/H1nI2MWbCJPYfSqFq6MD3PrUyPJpUoV8y64ozJU3/84da6GzUKihSBm26Ce++F6tW9riygLIwCLFzD6MDhdCYu2cLwORuZv3E3sdFRXFK/PL3PrUyLGqWJirJWkDGeWrQI3nzTDXbIyHCb+j3wQL4ZfWdhFGDhFEYZmcqva/5kzPxNTF66lUNpGdQsW4TezarQrXElTisS63WJxphjbdoE770HH30Eu3e7MHrgAbjiirCep2RhFGDhEEYrt+7j2wXJjF2wiW17D1OsYAyd48+ge+OKNKlayu4FGRMO9u+Hzz6Dt96CtWtdt13//nD99VC5co5vDzUWRgEWqmG0bW8q4xdtZsyCTSzbvJeYKKHdmWXp1rgSF9YtZ/sHGROuMjJg3Dg3LPynn9xGfxdeCH37Qrdu7j5TGLAwCrBQCqPt+w4zeekWJizeQsL6nahCfKUSdGtUkcsbnEHponFel2iMCaS1a+Hzz2HoUFi3DooWhR49XDC1aQNRobssl4VRgHkdRjv2H2by0q1MXLyF2et2kKlQu1xROsefQaf4CtQqV9Sz2owxeSQzE375xYXSyJFuEm21atC9O1x+udtbKSa0Fi22MAowL8Jo295Ufli+je+WbuW3tTvIyFRqlC1C5/gz6BxfgTqnF8vTeowxIeTgQRgzxrWYpk2DtDQ3V6lTJxdMHTtC8eJeV2lhFGh5EUaqSmLKfqYs38aU5dtYlLQbgGqlC9MpvgKd48+gbvliNhDBGPNPe/fClCnuHtOkSW61hwIFoF07F0zt20Pduu6+Ux6L6DASkY7AO7idXgep6svHvN4GtxNsPNBLVUfl9JnBCqOMTGVh0i6mLHMBtO5P9x0NKpXg4rPL0+Gs06ldrqgFkDHGP+np8NtvMH68C6eVK935smXd/aU2baBtWzjnnDy51xSxYSQi0cAqoAOQDCQAvVV1eZZrqgHFgQeBcXkdRjv2H2bm6u3MWLmdmau2s+tgGgWihRY1SrsAqne6LU5qjAmMNWtgxgyYOdONytuwwZ0vWRLOP9+FU7Nm0LBhULr1Qj2MgnmHrRmQqKprAURkONAV+CuMVHW977XMINbxl8xMZfGmPUz/I4UZq7azOHk3qlC6SCwX1C1HuzPL0e7MshQvWCAvyjHGRJKaNd3jppvc8YYNfwfTzJmuBXVU7drQqBE0buwejRpBmTLe1J1HghlGFYGkLMfJQPMgfl+2vknYyCvfrWTngSOIQMPKJbmvfR3anVmW+meUsOV4jDF5q2pVN4H2+uvd8datMH/+3485c2DEiL+vr1LFrTbeu7c39QZZaI09PAERGQAMAIiNPbUldMoVL0jbOmVpd2ZZ2tQuSylbiscYE0rKl4fLLnOPo3buhAULXDgtWOCuyaeCec+oJfCMql7iO34MQFVfOs61nwETvBzAYIwx+Vmo3zMK5hCOBKC2iFQXkVigFzAuiN9njDHmFIlIRxFZKSKJIvJoXn9/0MJIVdOBO4HvgRXACFVdJiLPiUgXABE5V0SSgauAj0VkWbDqMcYYc3y+0c/vA5cCZwG9ReSsPK3BJr0aY0z+l1033cncVgmW0F3VzxhjTCDFiMjcLI8BWV473ujninlaXF5+mTHGGM+kq2pTr4s4EWsZGWOM2QRk3TGwku9cnrEwMsYY4/noZ+umM8aYCKeq6SJydPRzNDBYVfN0dHPYjabzrWN36BTfHgOkB7CcUGA/U3iwnyk85LefKevPU0hVQ7Y3LOzCKDdEZG4o38A7FfYzhQf7mcJDfvuZwunnCdmUNMYYEzksjIwxxngu0sJooNcFBIH9TOHBfqbwkN9+prD5eSLqnpExxpjQFGktI2OMMSEoYsLI6+XRA01EBotIiogs9bqWQBGRyiIyXUSWi8gyEbnH65pyQ0QKisgcEVnk+3me9bqmQBGRaBFZICITvK4lEERkvYgsEZGFIjLX63oCQURKisgoEflDRFb4FkMNWRHRTedbHn0V0AG3AGAC0FtVl3taWC6ISBtgPzBMVet7XU8giEgFoIKqzheRYsA84Ipw/e8kIgIUUdX9IlIA+AW4R1V/97i0XBOR+4GmQHFV7ex1PbklIuuBpqr6p9e1BIqIDAV+VtVBvlUVCqvqbq/rOpFIaRk1AxJVda2qHgGGA109rilXVHUmsNPrOgJJVbeo6nzf8324fbDydOXgQFJnv++wgO8R9r/9iUgloBMwyOtazPGJSAmgDfApgKoeCeUggsgJI8+XRzcnR0SqAY2A2d5Wkju+7qyFQArwg6qG9c/j8zbwMJDpdSEBpMAUEZl3zNYK4ao6sB0Y4utOHSQiIbvlOEROGJkwIiJFgdHAvaq61+t6ckNVM1S1IW4V5GYiEtZdqiLSGUhR1Xle1xJg56lqY9xOp3f4usHDWQzQGPhQVRsBB4CQvlceKWHk+fLoxj++eyujgS9V9Vuv6wkUXxfJdKCj17XkUmugi+8ey3DgQhH5wtuSck9VN/n+TAHG4Lr2w1kykJylJT4KF04hK1LCyPPl0U3OfDf8PwVWqOqbXteTWyJSVkRK+p4Xwg2g+cPbqnJHVR9T1UqqWg33/9E0Vb3O47JyRUSK+AbM4OvKuhgI61GqqroVSBKRM32nLgJCeiBQRGwhEQrLoweaiHwNtAPKiEgy8LSqfuptVbnWGrgeWOK7zwLwuKpO8rCm3KgADPWN5owCRqhqvhgKnc+cDoxxvwsRA3ylqt95W1JA3AV86fsFfC3Qz+N6shURQ7uNMcaEtkjppjPGGBPCLIyMMcZ4zsLIGGOM5yyMjDHGeM7CyBhjjOcsjIwxxnjOwsgYY4znLIyMMcZ47v8Bi97FyQaZeEEAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "ax = plt.subplot(111)\n",
    "ax.plot(sol.t, sol.y.T[:, 0])\n",
    "ax2 = plt.twinx(ax)\n",
    "ax2.plot(sol.t, sol.y.T[:, 1], color='r')\n",
    "\n",
    "ax.set_ylabel('Biomass', color='b')\n",
    "ax2.set_ylabel('Glucose', color='r')"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.6"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": true,
   "title_cell": "Table of Contents",
   "title_sidebar": "Table of Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}