{
 "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": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Scaling...\n",
      " A: min|aij| =  1.000e+00  max|aij| =  1.000e+00  ratio =  1.000e+00\n",
      "Problem data seem to be well scaled\n"
     ]
    }
   ],
   "source": [
    "import cobra\n",
    "from cobra.io import load_model\n",
    "model = load_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": 3,
   "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:05, 34.89it/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]])\n",
       " y_events: [array([[0.87280538, 0.25178492]])]"
      ]
     },
     "execution_count": 7,
     "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": "iVBORw0KGgoAAAANSUhEUgAAAaMAAAD5CAYAAACK91rRAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6q0lEQVR4nO3dd3hU1dbH8e9KQkLvVULvgnSpCiigKCiKXIV7FSyI2MAulnttr+Xae0GKqAjSL01FBERRIfQqEHrognQSUtb7xx40IpCBzOTMZNbneeZh5kxbeRR+OfusvbeoKsYYY4yXorwuwBhjjLEwMsYY4zkLI2OMMZ6zMDLGGOM5CyNjjDGeszAyxhjjuZhgfriIdALeAqKBwar60knPVwKGAqWAfcCNqpp0ps+MiorSfPnyBaliY4zJnY4ePaqqGrInIBKseUYiEg2sBToCSUAC0FNVV2V6zRhgiqoOF5FLgVtU9aYzfW6BAgX0yJEjQanZGGNyKxE5qqoFvK7jdIKZks2ARFXdoKrHgVFA15Necz4w03d/1imeN8YYEwGCGUblga2ZHif5jmW2FOjmu38tUEhESpz8QSLSV0QWiMiCtLS0oBRrjDHGO16PHz4EtBWRxUBbYBuQfvKLVHWQqjZV1aYxMUG9zGWMMcYDwfyXfRtQIdPjeN+xP6jqdnxnRiJSELhOVfcHsSZjjDEhKJhnRglADRGpIiKxQA9gUuYXiEhJETlRw2O4zjpjjDERJmhhpKppwD3AN8BqYLSqrhSRZ0Xkat/L2gFrRGQtUAZ4Plj1GGOMOT0RGSoiu0VkRaZjxUXkWxFZ5/uzWNC+P9y2kLDWbmOMOXtZtXaLSBvgMPCpqtbzHXsZ2KeqL4nIQKCYqj4ajPoipxvgl1/gu++gQwdo2hSio72uyBhjyMhQ9h09zp5DKRxJSeNYajrHjqdzLDWd5D/uZ3AsNZ32tUvToELRoNShqnNEpPJJh7viRrAAhgOzAQujbPnhB3jySXcrWhQuuQQ6dnThVL06iHhdoTEmFzqYnMq6XYdI3H2YHQeS2X0ohd0HU9h9KJndB1P47XAKaRn+jVCVLhSXnTCKEZEFmR4PUtVBWbynjKru8N3fibucEhSRNUy3Zw/MnAkzZsC338Lmze54xYoumE7cihcPXMHGmIhwOCWNtbsOsW7XIdbuOuy7f5idB5P/8rriBWIpXSiO0oXzuj99t1KF8lIobwz5YqPJlyeavHmi/7ifL080cTFRREWd+y/N/qzA4DszmpJpmG6/qhbN9PzvqhqU60aRFUaZqcL69S6UZsxwIbV/P0RFQatWcOWV7la/vp01GWP+5sCxVBI27uOXDXv5ZeNeVm4/yIl/TvPmiaJ66YLULF2IGmUKUbNMQWqULkTZInmJjfFmeuc5htEaoJ2q7hCRcsBsVa0VlPoiNoxOlp4O8+fDV1/B1KmwaJE7Xr78n8HUvj0UKhT47zbGhLzThU9sTBRNKhajWZXi1CtfhJplChJfLD/R2TiLCYZzDKNXgL2ZGhiKq+ojQanPwug0duyAr7+GadNg+nQ4eBBiY901pmuugauvhjJBGz41xoSAA0dT+WbVTiYv3c5P6/eSnqF/hE+LqiVoUbU4DSoUJW+e0G+I8qObbiSuWaEksAt4CpgIjAYqApuB61V1X1DqszDyQ2oqzJ0LkyfDhAmwcaMbumvVCq691oVTtWo5W5MxJigOp6QxY9UuJi/dzpx1e0hNVyoWz0/n+uVoV7NU2ITPyUJ91W4Lo7OlCsuXw8SJLpiWLHHH69VzwXTDDVC3rnf1GWPOWmp6Bt+u2sWkJduZtWY3KWkZnFckL53rl+OqBudxQfkiSJhfO7YwCjDPw+hkmza5YJo40bWPZ2S4MLrhBnerWdPjAo0xp7P/6HG+mL+FT3/azM6DyZQsGEeX+uXoUr8cjSsWy1b3WqixMAqwkAujzHbuhHHj4MsvXTABNGwIPXrA9ddDlSqelmeMcRJ3H2bY3I2MW5REcmoGF1UvyS2tK9OuVumQazwIFAujAAvpMMosKQnGjHHBNG+eO9a8Odx0kwunEn/btskYE0Sqyo+JvzHkx43MXrOH2Jgorm1YnlsuqkztsoW9Li/oLIwCLGzCKLNNm2D0aBgxApYtgzx5oEsX6N0brrjCdekZY4JCVZm+ahevT1/Lml2HKFkwjl4tK/HP5hUpWTDO6/JyjIVRgIVlGGW2dCl8+qkLpl27oGRJ6NnTBVPjxjbB1pgAWrBpHy9+9SsLN/9OtVIFuLNdda5qUI64mPDrhssuC6MAC/swOiEtzc1fGj4c/vc/SElxjQ99+kCvXrYkkTHZkLj7EP/9eg3frtpF6UJx3N+xJv9oEk9MtNebW3vHwijAck0YZbZ/vxvGGzrUXV+Ki4N//AP69oWLLrKzJWP8tOtgMm/OWMuXCVvJHxtDv7ZVufWiKuSPjZw1oU/HwijAcmUYZbZ0KXz8MXz2mVv1oXZtF0q9elnTgzGncex4Ou/PTuTjHzaQnqH8q3kl7r20OiUi6JpQViyMAizXh9EJR464s6VBg9xeTHFx0L073Huv68ozxgDwU+JvDBy/nC37jnJVg/N46LKaVCoRsv/meiaiw0hEOgFvAdHAYFV96aTnK+I2bCrqe81AVZ12ps+MmDDKbNkyd7b06afubKlZMxgwwIWTdeKZCHXgWCovTlvNqIStVC6Rnxe71adlNRs9OJ2IDSMRiQbWAh2BJCAB6KmqqzK9ZhCwWFU/EJHzgWmqWvlMnxuRYXTCoUMukN5+G9auhbJl4c474Y47bNFWE1G+XrGT//xvBb8dTuH2NlW5v0PNsFwvLieFehgFs7WkGZCoqhtU9TgwCreFbWYKnJhtVgTYHsR6wl+hQnD33bB6tdvqolEjeOoptzlg795/bnthTC61+1Ayd41YSL/PF1KiYBz/u/siHruijgVRLhDMM6PuQCdV7eN7fBPQXFXvyfSacsB0oBhQAOigqgtP8Vl9gb4AsbGxTVJSUoJSc1haswbefReGDXPXmdq3h0cfdVtdWBeeySVUlbELk/i/qas5lprOgPY16NumKnkiuFX7bEXymZE/egKfqGo8cCXwmYj8rSZVHaSqTVW1aUyMtWj+Ra1a8M47sG0bvPKKO2u67DJo0sQtRZSW5nWFxmTLweRU7v5iEQ+PXUbNMgX5asDF3H1JdQuiXCaY/zW3ARUyPY73HcvsNtzGTajqz0Be3MZO5mwVKQIPPQQbNsDgwXD0qFsDr1Yt+OADOHbM6wqNOWtLt+6n89s/8M3KXQy8ojZf9m1JtVIFvS7LBEEwwygBqCEiVUQkFugBTDrpNVuA9gAiUgcXRnuCWFPuFxcHt90Gq1bB+PFQqhTcdRdUrgwvvOC68YwJcarKkB830v3Dn0hPV0bf0YJ+bavlqi0dzF8Fu7X7SuBNXNv2UFV9XkSeBRao6iRfB93HQEFcM8Mjqjr9TJ8Z0d1050IV5syB//7XNT0UKwYPPAD9+0Ph3L9SsQk/+48e56Exy5ixehcd6pTh1X/Up2h+m8KQXaF+zcgmvUaShQvhmWfc9unFisGDD7pJtBZKJkQs3Pw7/UcuZvehZAZeUYdbW1cO+x1WQ4WFUYBZGAXAwoXw9NMwZYpbkPVEKBUq5HVlJkKpKoPmbOCVb9ZQrmhe3u3ZmAYVinpdVq5iYRRgFkYBtGCBC6WpU10oPfSQG74rELL/v5pcKDk1nYfGLGXKsh1cUa8sL11XnyL58nhdVq5jYRRgFkZBkJDghu+mTnWrOjz9NNx6q9sE0Jgg2nUwmb6fLmDZtgM8cnlt+rWtasNyQRLqYWSN+gYuvNAN2f34I1SrBv36Qb16MHasa4AwJghWbDtA13fnsm73YT66sQl3tqtmQRTBLIzMn1q3hh9+cJv9xcS4PZVatIDZs72uzOQy05bvoPuHPxEdJYzt14rL6pb1uiTjMQsj81cicPXVbl+lIUPcyg6XXAJXXulWDzcmG1SVt79bx10jFnF+ucJMvLs1559n3ZzGwsicTkyMu260bp2bo/Tzz9Cwodvob4/NSzZnLzk1nf6jlvD6t2u5tlF5vri9BaUK2eZ3xrEGBuOf33+HZ591i7IWKOBWC7/7bttPyfjlt8Mp3DZ8AUu37ufhy2txl10fynGh3sBgYWTOzurVcP/98M03bt27N9+ETp28rsqEsK37jnLTkHnsPJjMmzc0olM9uz7khVAPIxumM2enTh23rNCUKZCRAVdcAV26uM3+jDnJ6h0H6fbBT/x+NJURfZpbEJnTsjAyZ08EOneGFSvg1Vfd2nf16rlJs4cOeV2dCRHzNuzl+o9+JlqEMf1a0qRSca9LMiHMwsicu9hYt5TQunXQqxe8/jqcfz5MmGDzkyLc9JU7uWnofEoVimPcXa2oWcaWmjJnZmFksq9MGbeH0k8/uWWFunWDrl1h82avKzMeGJ2wlX6fL6ROucKM7deK8kXzeV2SCQMWRiZwWrRwi7C++ip89507S3rlFUhN9boykwNUlfdnJ/LIuGVcVKMUX/RpTvEC1m0ZLkTkfhFZKSIrRGSkiOTNye+3MDKBFRPjhu5Wr4YOHeCRR9wW6D//7HVlJogyMpT/m7qal79ew9UNzmNwr6YUiIvxuizjJxEpD/QHmqpqPdwedD1ysgYLIxMcFSu6ZYUmTHBzlFq3dmveHTjgdWUmwDIylMcnLGfIjxu5uVVl3ryhIbEx9k9LGIoB8olIDJAf2J6TXx7U/2NEpJOIrBGRRBEZeIrn3xCRJb7bWhHZH8x6jAeuuebPuUkff+y67r7+2uuqTICkZygPjV3KqISt3HtpdZ666nzbGjx0xYjIgky3vieeUNVtwKvAFmAHcCCrXbcDLWiTXkUkGlgLdASSgASgp6quOs3r7wUaqeqtZ/pcm/QaxubPh1tugVWr3J+vvw5Fi3pdlTlHaekZ3D96KZOXbueBjjXp376G1yWZMzjTpFcRKQaMA24A9gNjgLGq+nlO1RfMM6NmQKKqblDV48AooOsZXt8TGBnEeozXmjWDRYvg8cfh00+hbl23h5IJO8fTMrh35GImL93OwCtqWxCFvw7ARlXdo6qpwHigVU4WEMwwKg9szfQ4yXfsb0SkElAFmHma5/ueOLVMS0sLeKEmB8XFwfPPw7x5rg28Sxe4+WZ3XcmEhZS0dO4asZCvVuzk313Op1/bal6XZLJvC9BCRPKLWzSwPbA6JwsIlauMPXCnhOmnelJVB6lqU1VtGhNjHTq5QpMmbtvzJ5+Ezz93Z0mTJ3tdlclCcmo6fT9dyIzVu3mua11uu6iK1yWZAFDVecBYYBGwHJcNg3KyhmCG0TagQqbH8b5jp9IDG6KLPHFx8Nxz7lpSqVJuH6U+feDwYa8rM6dw7Hg6fYYvYM66PbzU7QJualnZ65JMAKnqU6paW1XrqepNqpqSk98fzDBKAGqISBURicUFzqSTXyQitYFigE1EiVSNG0NCAjz2GAwdCo0auWE8EzKOpKRx87D5/LT+N17t3oAezSp6XZLJZYIWRqqaBtwDfIMbexytqitF5FkRuTrTS3sAozTc9rIwgRUbCy+84LY4P37czUt67jmwa4SeO3Y8nduGJ5CwaR9v3NCQ65rEe12SyYVsPyMTevbvdxv3ffEFtGrlrilVsWsTXkhOdUNzP63/jTduaEjXhqfsQTJhwPYzMuZsFS0KI0a424oV0KCBawUPs1+cwl1yajp3fLaQuet/45XuDSyITFBZGJnQ9c9/wrJl0LAh9O4NPXq4syYTdMfTMrhrxCK+X+uaFWxozgSbhZEJbZUqwaxZ7nrS+PGuJXzhQq+rytVS0zO454tFzPx1N/93TT1uuNCaFUzwWRiZ0Bcd7Trt5sxx21G0agXvvmvDdkGQlp7BgFGLmb5qF89cXZcbW1TyuiQTISyMTPho2RIWL4aOHeHee+H6620V8ABKz1AeGL2Uact38mTnOvRuVdnrkkwEsTAy4aVECZg0CV5+2W1P0bixW+/OZEt6hvLwmKVM8q011+fiql6XZCKMhZEJP1FR8PDDbtju+HF3xvT++zZsd45UlScnLmf84m082LGmrTVnPGFhZMJXq1awZInbUfbuu+GGG+DQIa+rCiuqbofWkfO3cvcl1bjXVt82HrEwMuGtRAm3wOpLL7luu+bNYc0ar6sKG2/MWPfHDq0PXVbL63JMBLMwMuEvKgoefRS+/Rb27HH7JtkK4Fn66Pv1vP3dOq5vGs9/upyP2znAGG9YGJnc45JL3BykGjXcCuDPPAMZGV5XFZI++2UzL371K13ql+PFbvVtq3DjOQsjk7tUrAg//AC9esHTT8M111j790nGLUzi3xNX0KFOad64oSHRFkQmBFgYmdwnXz745BN45x346is3bLdqlddVhYSvlu/g4bFLaV29BO/+szF5ou2fABMa7P9EkzuJwD33wMyZbj275s1dg0MEm7VmN/1HLaZRxWIMuqkpefNEe12SMX+wMDK528UXu+tI558P110Hzz4bkfOR5m3YS7/PFlKzTCGG3nwhBeJivC7JmL+wMDK5X3y8myDbuzc89RT8619w7JjXVeWY5UkHuG34AuKL5ePTW5tRJF8er0sy5m8sjExkiIuDYcPgxRdh5EjXebdzp9dVBd26XYfoNXQeRfPnYUSfFpQoGOd1ScacUlDDSEQ6icgaEUkUkYGnec31IrJKRFaKyBfBrMdEOBEYONBdO1q+3DU2LF3qdVVBs3XfUW4cMo+Y6Cg+v605ZYvk9bokY04raGEkItHAe8AVwPlATxE5/6TX1AAeA1qral3gvmDVY8wfrr0WfvzRzUFq3TpXTpDdfTCZG4fMIzk1g89ua0blkiG727QxQHDPjJoBiaq6QVWPA6OArie95nbgPVX9HUBVdwexHmP+1KgRzJ8PdepA167w6qu5prFh/9Hj3DRkPnsOpTDslgupXbaw1yUZk6VghlF5YGumx0m+Y5nVBGqKyFwR+UVEOp3qg0Skr4gsEJEFaWlpQSrXRJzzzoPvv4fu3d0q4H36uFXAw9jhlDR6D0tg429H+LhXUxpXLOZ1Scb4xev+zhigBtAOiAfmiMgFqro/84tUdRAwCKBAgQK549dXExry54dRo6B2bXjuOdiyBcaNg8LhdzaRnJpO308XsGLbAd7/V2NaVy/pdUnG+C2YZ0bbgAqZHsf7jmWWBExS1VRV3QisxYWTMTknKsrNPxo2DGbPdnOTtp38v2poS03P4N6Ri/lp/V5e6V6fy+uW9bokY85KMMMoAaghIlVEJBboAUw66TUTcWdFiEhJ3LDdhiDWZMzp3XwzTJkCGza4DftWrvS6Ir9kZCiPjF3Gt6t28czVdenWON7rkow5a0ELI1VNA+4BvgFWA6NVdaWIPCsiV/te9g2wV0RWAbOAh1V1b7BqMiZLl1/uJsimprpOu9mzva7ojFSVpyevZIJvl9berSp7XZIx50Q0zDqIChQooEeOHPG6DJPbbd4MV1wB69fD8OHQo4fXFZ3Sa9PX8M7MRG6/uAqPX1nH9iQypyUiR1U1ZHv8bQUGY06lUiWYOxdatICePUOy9fvjORt4Z2YiNzStYEFkwp6FkTGnU6wYfPMNXH+9a/0eMADS072uCoBR87fw/LTVdL6gHC90u8CCyIQ9r1u7jQltefO6tewqVIDXXoNdu+CzzyA21rOSpi7bwWMTltO2ZinbHM8ElIgUBQYD9QAFblXVn3Piuy2MjMlKVJQbpitb1p0h7d/v5iIVLJjjpcxes5v7vlxMk4rF+PDGJsTG2OCGCai3gK9VtbuvCzp/Tn2xNTAYczaGDXMrNVx4IUydCiVK5NhXJ2zax01D5lG1ZEFG9m1hW0GYs5JVA4OIFAGWAFXVg2Dw69cqEaqJEOe7306E/iIUDWplxoSiW25xZ0VLlkCbNjk2OXbFtgPc+kkC5xXJx6e32Z5EJiiqAHuAYSKyWEQGi0iOdd/5e44/DkgXoTpuWZ4KgG33YCLTNdfAV1/B1q1uLtK6dUH9usTdh+g1dD6F8+bhsz7NKWl7EplzE3NijU/fre/JzwONgQ9UtRFwBDjl1j9nJFIJkQ6++/kQKeTP2/wNowxV0oBrgXdUeRgod9ZFGpNbXHIJzJoFR464QFq8OChfs3XfUW4cPJ8oET7v05zyRfMF5XtMREhT1aaZboNOej4JSFLVeb7HY3Hh5D+R233v+8h3JB630k6W/A2jVBF6Ar2BKb5jNk5gIluTJm5fpHz5oG1btwJ4AO06mMy/Bs/jWGo6n/dpRhXbk8gEkaruBLaKSC3fofbAqrP8mLuB1sBB34euA0r780Z/w+gWoCXwvCobRagCfHaWRRqT+9Sq5SbHxse7pYSmTMn6PX7Yd+Q4Nw6ex97DKQy/tZntSWRyyr3ACBFZBjQEXjjL96fg9q9zRGJwLeJZOutuOhGKARVUWXZWbwwQ66YzIWnvXujUyTU2jBzp9kg6RweTU/nXx/NYu+sQn9zSjJbVcq5jz+ReObIckMjLwH6gFy7Y7gJWofpEVm/1t5tutgiFRSgOLAI+FuH1c6/YmFymRAmYMQOaN4cbbnATY8/BsePp3PZJAqt3HOTDG5tYEJlwMxDXkbccuAOYBjzpzxv9HaYrospBoBvwqSrNgQ7nUKgxuVeRIm75oHbtoHdv+OijLN+SWUpaOn0/W8DCzb/zVo9GXFLbr6F2Y0KHagaqH6P6D6AvMA8/h9/8DaMYEcoB1/NnA4Mx5mQFCrjrRldcAf36wRtv+PW2tPQMBoxcwg/rfuOlbvXpXN+aVU0YEpmNSGFEigMLgY8R8esvgb9h9Cxu76FEVRJEqAoEd3KFMeEqXz6YMAGuuw4eeACef/6ML0/PUB4eu4yvV+7kP13O5/oLK5zx9caEsCKo/jGKhmpzXFdelvwKI1XGqFJflbt8jzeoct05l2tMbhcbC6NGwY03wpNPwuOPn3ILiowM5YkJy5mweBsPX16LWy+q4kGxxgRMDCLnNIrm10KpIuQFbgPqAnlPHFfl1rP5MmMiSkyM25gvXz548UU3QfbNN8G33YOq8szklYxK2Er/S6tz9yXVva3XmOw7MYo2F9UERPweRfN3mO4zoCxwOfA9blbtoazeJCKdRGSNiCSKyN+WlRCRm0Vkj4gs8d36+FmPMeEhKso1MgwYAG+/DXffDRkZqCovffUrw3/ezO0XV+H+jjW9rtSY7FMdg2p9VO/0Pd6Aql+jaP5uIVFdlX+I0FWV4SJ8AfxwpjeISDTwHtARt8xEgohMUtWTZ/R+qar3+FmHMeFHxDUy5M0L//0vpKfzRrf7+WjOBnq1rGS7tJrcQyQeeAe3CgO4nBiAalJWb/V7OSDfn/tFqAcUIeslHpoBiaq6Qd2M3FFAVz+/z5jcRcQN1T3+OAwaRLlHBtCjcXmevqquBZHJTYYBk4DzfLfJvmNZ8jeMBvlWXvi374tWAS9n8Z7ywNZMj5N8x052nYgsE5GxInLKNiIR6Xtipdm0tDQ/SzYmxIgwpNNtvNWqBz2XTeeFr94iSjO8rsqYQCqF6jBU03y3T4BS/rzR3266war8rsr3qlRVpbQqH2anYp/JQGVVrQ98Cww/9ffroBMrzcbE2Oa0JjyNmLeZ56au5tc7HyL9P08R9ckncOutkJ7udWnGBMpeRG5EJNp3uxHY688b/e2mK4pba6hy5veo0v8Mb9uG2/fohHjfsT+oauYiB5P12ZYxYWn0gq08MWEF7WuX5q0ejYiOaQIx0fCf/0BGBnzyCURHe12mMdl1K+6a0Ru4BVJ/wi20nSV/TzOmAb/g1hvyd1whAaghIlVwIdQD+GfmF4hIOVXd4Xt4NbDaz882JmyMWbCVR8ct4+IaJXnvX42JjfENSPz73y6AnnjCBdLw4a4d3JhwpboZ92/5WfP3//y8qjxwNh+sqmkicg+u5zwaGKqqK0XkWWCBqk4C+ovI1UAasA+4+Wy+w5hQN25hEo+MW8ZF1Uvyca+m5M1z0tnP44+7QBo40AXSZ59ZIJnwJTIc1z233/e4GPAaqlnOSfVrCwkR7gcO42bUppw4rsq+c6v43NkWEiZcTFicxAOjl9KqWgmG9L7w70GU2auvwsMPQ48eFkgmKHJoC4nFuC3Lz3zsFPz9P/448ArwBH9ulKRA1bMo05iI8b8l23hw9FJaVCnB4F5ZBBHAQw+55YIeecRNlP30U7uGZMJRFCLFUP0dwLdgql85428YPYib+PrbudVnTOSYtHQ793+5hGZVijPk5qbki/UzVB5+2HXWPfaYC6JhwyyQTLh5DfgZkTGAAN2BM68U7ONvGCUCR8+tNmMix+Sl27lv1GKaVi7O0JsvJH/sWQ63DRzoAunJJ10QDRnizpSMCQeqnyKyALjUd6Qbf19155T8/ZtyBFgiwiz+es3oTK3dxkSUqct2cN+XS2haqTjDziWITnjiCUhLg6efdoE0aJAFkgkPIhVx/QWT/nJMdUtWb/X3b8tE380YcwpTlm1nwKglNKpQlKG3XEiBuGw2IDz1lDtDeu45F0gffGCBZMLBVP7sK8gHVAHW4HZ8OCO//sb4FkeNBU4sLbxG9Y/16oyJaBMWJ/Hg6KU0qVSMYbc0o2B2g+iEZ55xgfTCCy6Q3nvvj+0njAlJqhf85bFIY3D74GXF3xUY2uGW6tmEuyhVQYTeqsw5mzqNyW1GJ2zl0fHLXNdc76bZPyPKTAT+7//ckN3LL7szo3fesUAy4UN1ESLN/Xmpv39zXgMuU2UNgAg1gZFAk3Or0Jjw9/kvm3ly4gourlGSQTedRdfc2RCBl15yZ0ivveZ2kH3tNQskE5pEMi+OEAU0Brb781Z/wyjPiSACUGWtCHn8r9CY3GXY3I08M3kVl9Yuzfv/apz1PKLsEIFXXnFnSG+8AXFxbujOAsmEnkKZ7qfhriGN8+eN/obRAhEGA5/7Hv8LWOB3ecbkIoPmrOeFab9yed0yvNMz01pzwXRig76UFHemFBfnuu2MCSWqz5zrW/0NozuBu+GPVu4fgPfP9UuNCVfvzlzHq9PX0qV+Od64oSF5onOww03ENTEcP+6aG2Jj3dp2xnhNZDJ/dtH9nWqWi6f6202XArzuuxkTcVSVN2as4+3v1tGtUXle7l6fmJwMohOioty8o+PH3XykuDh48MGcr8OYv3r1FMdOhJNf48lnDCMRRqtyvQjLOUXqqVLfny8xJpypKs9PXc3gHzdyfdN4XuxWn+goD6/XnFgqKCXFrWkXGwv33utdPcZAUSAe1fcAEJmP2+FVgUf9+YCszowG+P7scm71GRPe0tIzeHzCckYvSKJ3y0o8dVVdorwMohNiYmDECEhNhf793RlS375eV2Ui1yO4PetOiAWaAgWAYcCYrD7gjGGkyg7fn5tPHBOhJLBX9Qzjg8bkAilp6QwYuYSvV+6kf/sa3N+hBhJKHWx58sCoUdCtG9xxhztDuvlmr6sykSkW1a2ZHv+I28l7LyJ+bVtxxkFvEVqIMFuE8SI0EmEFsALYJUKnc6/bmNB29HgafYYv4OuVO/l3l/N5oGPN0AqiE+LiYNw46NgRbr0VvvjC64pMZCr2l0eq92R6VMqfD8jqCuy7wAu4Ca4zgT6qlAXaAC9m9eEi0klE1ohIoogMPMPrrhMRFZGm/hRtTDAdOJrKjYPnMTfxN17pXp/bLqridUlnljcvTJwIbdtCr14unIzJWfMQuf1vR0XuAOb78wFn3OlVhCWqNPTdX61KnUzPLVbltLv3iUg0sBboCCQBCUBPPWk5cREphJsYFQvco6pnnL9kO72aYNp9KJleQ+azYc8R3u7ZkE71ynldkv8OH4bLL4f582H8eLjqKq8rMiHEn51eff9uLwC2qar/vQIipXGLaacAi3xHmwBxwDWo7srqI7I6M8rIdP/YSc9ldc2oGZCoqhtU9TgwCuh6itc9B/wXSM7i84wJqq37jvKPD39my76jDL35wvAKIoCCBWHaNGjYELp3h+nTva7IhJ8BwOqzfpfqblRb4f493+S7PYtqS3+CCLIOowYiHBThEFDfd//E4wuyeG95IPMFrSTfsT+IW9G1gqpOPdMHiUhfEVkgIgvS0tKy+Fpjzt7aXYfo/uFP7D+ayud9mnNRjZJel3RuihSBb76BOnWga1eYPdvrikyYEJF4oDMw+Jw/RHUmqu/4bjPP5q1nDCNVolUprEohVWJ89088ztbadCIShZtEm+WMPVUdpKpNVbVpTEwAV0U2Bvhlw166f/ATqvDlHS1oXLFY1m8KZcWLw7ffQtWq0KULzJ3rdUUmNMSc+KXedzt5LsCbuBbtjL+/NfiCOYV8G1Ah0+N437ETCgH1gNkisgloAUyyJgaTk6Ys206vIfMpXTgv4+9qRe2yhb0uKTBKlYLvvoPy5eGKK9x1JBPp0k78Uu+7DTrxhIh0AXar6kKvigtmGCUANUSkiojE4iZE/bEVraoeUNWSqlpZVSsDvwBXZ9XAYEygDP1xI/eOXEz9+CKM7deS+GL5vS4psMqWdYFUsqRrbFiyxOuKTOhqDVztOzEYBVwqIp+f+S2BFbQwUtU04B7gG9wFsdGqulJEnhWRLBfNMyZYMjKUF6at5tkpq7j8/LJ83qc5RfPHel1WcMTHw8yZUKgQdOgAK1Z4XZEJQar6mKrG+04MegAzVfXGnKzhjK3dochau012pKSl8/CYZUxaup1evuV9PF1nLqckJkKbNm6Tvu+/h9q1va7I5DB/Wrt9r2sHPHRWrd0BYGFkIsbB5FT6fbaQn9bv5dFOtenXtmporqoQLL/+Cu3auZW/v/8eatTwuiKTg/wNI694sAa+MTlvx4FjXP/hz8zfuI/Xr2/Ane2qRVYQgTsbmjHDLa566aWwcaPXFRnzBwsjk+stTzrANe/NZeu+owy75UK6NY73uiTv1KvnAunIEbjkEtiyxeuKjAEsjEwu9/WKHfzjo5+IiYpi3F2tuLiGX2s25m4NGrh5SPv3uzOkbduyfIsxwWZhZHIlVeX92Yn0+3wRdcoVZuLdrXPPHKJAaNLErdSwe7cLpJ07va7IRDgLI5PrHE/L4OGxy3j56zVc1eA8Rt7eglKF4rwuK/Q0b+7Wstu2Ddq3hz17vK7IRDALI5Or/H7kODcOmcfYhUkMaF+Dt3s0JG+eaK/LCl0XXQRTprhmhg4dYO9erysyEcpau02usX7PYW79JIEdB5J5pXt9ujYsn/WbjDNjhlvHrk4dd79ECa8rMgFmrd3G5IA5a/dw7XtzOZycxsjbm1sQna0OHWDSJFi92u0au2+f1xWZCGNhZMKaqvLh9+u5edh8yhXJx8S7W9OkUnGvywpPl10G//sfrFrlwskCyeQgG6YzYevo8TQeHruMqct20Ll+OV6+rj4F4myLkWz7+mu3F1K9eq4FvLiFe24Q6sN0FkYmLG3ee4Q7PlvI2l2HeKRTbe5oE2FL+wTbV1/BNdf8OUm2WJjv8WQsjALNwsh8v3YP/UcuBuCdno1oU9MmsgbFtGlw7bVwwQUukIoW9boikw2hHkZ2zciEDVXlg9nruWXYfMoVycvkey6yIAqmK6+E8eNh+XLX1LB/v9cVmVzMzoxMWDiSksYjY5cxdbm7PvRK9/rkj7XrQzliyhTo1g0aNoTp0+0MKUyF+pmRhZEJeb/uPMhdIxax6bcjPNqpNn3t+lDOmzwZrrvODdlNn27zkMJQqIeRDdOZkKWqjE7YyjXvzeVQchqf92nOHW0jcOuHUHDVVTBxIqxc6Vb73r3b64pMLhPUMBKRTiKyRkQSRWTgKZ7vJyLLRWSJiPwoIucHsx4TPo4eT+PBMUt5ZNwymlQqxrT+F9OqWkmvy4psV17phuwSE90mfTt2eF2RyUWCNkwnItHAWqAjkAQkAD1VdVWm1xRW1YO++1cDd6lqpzN9rg3T5X5rdx3irhGLWL/nMAPa1+DeS2tExtbg4WLOHOjcGcqWhZkzoUIFrysyfojkYbpmQKKqblDV48AooGvmF5wIIp8CQHhdwDIBN3ZhEl3fncv+o8f5/Lbm3NehpgVRqGnTxl032r3b3bcdY00ABDOMygNbMz1O8h37CxG5W0TWAy8D/U/1QSLSV0QWiMiCtLS0oBRrvHX0eBoPj1nKQ2OW0qBCEab1v5jW1W1YLmS1bAnffQcHDkDbtrBundcVmTDneQODqr6nqtWAR4EnT/OaQaraVFWbxsRYO29uszzpAF3e/pGxi5Lof2l1RvRpQenCeb0uy2SlaVM3THfsmAuk1au9rsiEsWCG0TYg82ByvO/Y6YwCrgliPSbEpGco781K5Nr353IsNZ0RfZrzwGW1bFgunDRsCLNnQ0aGC6SlS72uyISpYIZRAlBDRKqISCzQA5iU+QUiUiPTw86AnetHiKTfj9Lz41945Zs1XF6vLF8PaGPdcuGqbl3X1BAX5wLphx+8rsiEoaCNealqmojcA3wDRANDVXWliDwLLFDVScA9ItIBSAV+B3oHqx4TOv63ZBtPTlyBKrz2jwZ0a1ze5g6Fu5o1Ye5ctw3FZZfBmDFusz5j/GQrMJgcc+BYKv+euIJJS7fTtFIx3rihIRWK5/e6LBNIe/a4+UiLF8PQodCrl9cVGZ9Qb+22bgCTI35K/I2Hxy5j58FkHuxYkzvbVSMm2vP+GRNopUq5poZrroHevd0Gfffd53VVJgxYGJmgOpScyotf/coX87ZQpWQBxvZrSaOKtjdOrlaokNt+4p//hPvvh99+g+eeAxuKNWdgYWSCZs7aPTw2fjnbDxzj9our8EDHWuSLjfa6LJMT4uJg9Gi48054/nkXSO+9B9H239+cmoWRCbiDyak8P2U1Xy7YStVSBRjbrxVNKtnZUMSJjoaPPoKSJeHFF2HvXvj8cxdUxpzEwsgE1Kw1u3l8/HJ2HUzmjrZVub9DTfLmsd+GI5YIvPCCC6QHH4Rdu2DCBNuCwvyNddOZgNh/9Dj/N3U1YxcmUaN0QV7uXt+uDZm/GjUKbr4ZKlaEqVOhRo0s32ICJ9S76SyMTLaoKhMWb+P5qavZfyyVO9pUpX/7GnY2ZE5t7lzo2hVU3f5IF1/sdUURw8IowCyMQsf6PYd5csIKft6wlwYVivLCtfWoe14Rr8syoS4x0W1BsWkTDBvmuu5M0GUVRiJSAfgUKIPbQWGQqr6VY/VZGJmzlZyazvuzEvnw+w3E5Yni0U616dmsoq0pZ/y3bx9ce61bRujZZ+HJJ631O8j8CKNyQDlVXSQihYCFwDWZ96ALJmtgMGflh3V7+PfEFWzae5SuDc/jic51KF3IVtg2Z6l4cbcn0u23w3/+486WPv4YYmO9rixiqeoOYIfv/iERWY3b9sfCyISOXQeTeX7qaiYt3U6VkgX4/LbmXFTDFjY12RAXB8OHQ/Xq8NRTsHkzjBtnnXbBEyMiCzI9HqSqg071QhGpDDQC5uVEYWDDdCYLyanpDPlxI+/NSiQtXbmzXTXubFfNGhRMYI0YAbfeCuXKwfjx0Lix1xXlOv42MIhIQeB74HlVHR/8ynzfa2FkTkVV+XrFTp6ftpqk349xed0yPH5lHSqVCNlmHBPu5s+H665zi61+8AHccovXFeUq/oSRiOQBpgDfqOrrOVOZ77stjMzJVm0/yLNTVvLLhn3UKlOI/1x1vm0BbnLGnj3Qo4dbbPWOO+Ctt2zFhgDxo4FBgOHAPlW9L8cKO/H9FkbmhL2HU3jt27WMmr+FIvny8MBlteh5YQVbXdvkrLQ0eOIJePllaN4cxo6F+Hivqwp7foTRRcAPwHIgw3f4cVWdliP1WRiZ5NR0PvlpE+/NSuTY8XRualmJ+9rXpEj+PF6XZiLZuHFuxYZ8+eDLL+GSS7yuKKxF9KRXEekEvIXb6XWwqr500vMPAH2ANGAPcKuqbj7TZ1oYBU5aegZjFybx5ox17DyYzCW1SvFE5zpUL13I69KMcVavdvOREhPhpZfc+nY2H+mcRGwYiUg0sBboCCQBCUDPzBOoROQSYJ6qHhWRO4F2qnrDmT7Xwij7VJVvVu7ilW9+Zf2eIzSqWJSBnWrTvKq11JoQdPCga2YYP96t3DBkCJQp43VVYSfUwyiYFwOaAYmqukFVjwOjgK6ZX6Cqs1T1qO/hL4ANDAfZLxv20u2Dn+j3+UIAPryxCePvbGVBZEJX4cLuutFbb8GMGXDBBTB5stdVmQALZhiVB7ZmepzkO3Y6twFfneoJEekrIgtEZEFaWloAS4wcK7Yd4JZh8+kx6Bd27E/mv9ddwDf3taFTvbKIDXuYUCcC/fvDwoVw3nlw9dXQrx/YKEmuERIrMIjIjUBToO2pnvfNEh4EbpguB0sLe8uTDvDWd2uZsXo3hfPGMPCK2tzcqrJNWjXhqW5dmDcP/v1vePVVmDXLTZht2tTrykw2BTOMtgEVMj2O9x37CxHpADwBtFXVlCDWE1GWJe3nrRnr+O7X3RTJl4cHO9akd+vKFM5rHXImzMXFubbvK66AXr2gZUt4+mkYONC2NQ9jwWxgiME1MLTHhVAC8E9VXZnpNY2AsUAnVV3nz+daA8OZLd26n7e+W8dMXwjdfnEVereqTCELIZMb/f473HWX27ivdWsYOhRq1vS6qpAU6g0MwW7tvhJ4E9faPVRVnxeRZ4EFqjpJRGYAF+BbKRbYoqpXn+kzLYxObeHm33lvViIzf91N0fx5uP3iqvRqWclCyESGL75woXTsGDz2mDtLymuryWcW0WEUDBZGf8rIUGb+upuP5qwnYdPvf4RQ71aVKRgXEpcDjck5O3e6eUhffOG2NH//fejQweuqQoaFUYBZGEFKWjr/W7KdQXM2kLj7MOWL5uP2i6tw/YUVyB9rIWQi3IwZ7ixp3Tq3i+xrr0HZsl5X5TkLowCL5DA6lJzKF/O2MHTuRnYdTKFOucL0a1uVzheUs/XjjMksOdmt2PDii245oRdfhL59I7rBwcIowCIxjLbsPcqnP2/iy4StHEpJo3X1EtzRphoX1yhpc4SMOZO1a91Z0nffQbNm8OabrvsuAlkYBVikhJGq8mPibwz/aRPf/bqbaBE61StL3zZVqR9f1OvyjAkfqjBypLuetHMndO0Kzz/v5ixFEAujAMvtYXQkJY3xi5IY/vNmEncfpkSBWP7ZvCL/al6JskWsO8iYc3bkiFtS6L//hUOH3BylZ56BSpW8rixHWBgFWG4No/V7DjPily2MWbiVQ8lp1I8vws2tKtO5fjniYiJ3nNuYgNu7111Peucdd9Z0553w+ONQurTXlQWVhVGA5aYwSk5N5+sVO/li/hbmb9xHTJRw5QXluLl1ZRpVKGrXg4wJpqQkd2Y0dCjkz++G8QYMgGLFvK4sKCyMAiw3hNGanYcYOX8LExZv48CxVCqVyM8NF1age5N4SheyoThjctSvv7q17saOhQIF4Lbb4L77oEoVrysLKAujAAvXMDqSksbU5TsYNX8Li7bsJzY6isvrlaXnhRVoUbUEUVF2FmSMp5Yuhddfd80O6eluU78HH8w13XcWRgEWTmGUnqH8tP43JizaxlcrdnIsNZ1qpQrQs1lFujWOp3iBWK9LNMacbNs2ePdd+PBD2L/fhdGDD8I114T1PCULowALhzBas/MQ4xcnMXHxNnYdTKFQ3hi61D+P6xqXp0mlYnYtyJhwcPgwfPIJvPEGbNjghu369IGbboIKFbJ8e6ixMAqwUA2jXQeTmbx0OxMWb2Pl9oPERAntapWiW+N4Lq1d2vYPMiZcpafDpEmuLfz7791Gf5deCr17Q7du7jpTGLAwCrBQCqM9h1L4asUOpizbQcKmfahC/fgidGtUnqsanEeJgnFel2iMCaQNG+Czz2D4cNi4EQoWhO7dXTC1aQNRobssl4VRgHkdRnsPp/DVip1MXbaDeRv3kqFQo3RButQ/j871y1G9dEHPajPG5JCMDPjxRxdKY8a4SbSVK8N118FVV7m9lWJCa9FiC6MA8yKMdh1M5ttVu/h6xU5+3rCX9AylaqkCdKl/Hl3ql6NmmUI5Wo8xJoQcPQoTJrgzppkzITXVzVXq3NkFU6dOULiw11VaGAVaToSRqpK4+zDTV+1i+qpdLN26H4DKJfLTuX45utQ/j9plC1kjgjHmrw4ehOnT3TWmadPcag958kC7di6YOnSA2rXddaccZmEUYMEKo/QMZcnW35m+0gXQxt/cdzSIL8JldcvS8fwy1Chd0ALIGOOftDT4+WeYPNmF05o17nipUu76Ups20LYtXHBBjlxriugwEpFOwFu4bccHq+pLJz3fBrcteX2gh6qOzeozAxlGew+nMGfdHmav2cOctXv4/WgqeaKFFlVLuACqU8YWJzXGBMb69TB7NsyZ47ryNm92x4sWhYsvduHUrBk0bBiUYb2IDSMRiQbWAh2BJCAB6KmqqzK9pjJQGHgImBTsMMrIUJZtO8CsX3cze+0eliXtRxVKFIilba1StKtVmna1SlE4b55z+nxjjPHb5s1/BtOcOW5n2hNq1IBGjaBxY3dr1AhKlszW10VyGLUEnlbVy32PHwNQ1RdP8dpPgCnBDKMvE7bw36/XsO/IcUSgYYWiXOILn3rnFbHleIwx3tq5ExYt+vO2eDFs2vTn8xUrutXGe/Y8p48P9TAKZu9heWBrpsdJQPNz+SAR6Qv0BYiNPbcldEoXzkvbmqVoV6sUbWqUopgtxWOMCSVly8KVV7rbCfv2uVA6EU5ly3pXX5CFViP8aajqIGAQuDOjc/mMS2qV5pJauXu/EmNMLlO8OLRv7265XDBbOLYBmRdwivcdM8YYE2JEpJOIrBGRRBEZmNPfH8wwSgBqiEgVEYkFegCTgvh9xhhjzoGv4ew94ArgfKCniJyfkzUELYxUNQ24B/gGWA2MVtWVIvKsiFwNICIXikgS8A/gIxFZGax6jDHGnFYzIFFVN6jqcWAU0DUnCwjqNSNVnQZMO+nYfzLdT8AN3xljjAmuGBFZkOnxIN/1eAhgw9m5CosGBmOMMdmWpqpNvS7idEJ3vXNjjDE5xfOGMwsjY4wxnjec2TCdMcZEOFVNE5ETDWfRwFBVzdGGsrBbtVtEMoBj5/j2GCAtgOWEAvuZwoP9TOEht/1MmX+efKoasqNhYRdG2SEiC0L5At65sJ8pPNjPFB5y288UTj9PyKakMcaYyGFhZIwxxnORFkaDsn5J2LGfKTzYzxQectvPFDY/T0RdMzLGGBOaIu3MyBhjTAiyMDLGGOO5iAkjr/fqCDQRGSoiu0Vkhde1BIqIVBCRWSKySkRWisgAr2vKDhHJKyLzRWSp7+d5xuuaAkVEokVksYhM8bqWQBCRTSKyXESWnLSYaNgSkaIiMlZEfhWR1SLS0uuaziQirhn59upYC3TErUabAPRU1VWeFpYNItIGOAx8qqr1vK4nEESkHFBOVReJSCFgIXBNuP53EhEBCqjqYRHJA/wIDFDVXzwuLdtE5AGgKVBYVbt4XU92icgmoKmq/uZ1LYEiIsOBH1R1sG+Jn/yqut/jsk4rUs6MPN+rI9BUdQ6wz+s6AklVd6jqIt/9Q7h9sMp7W9W5U+ew72Ee3y3sf/sTkXigMzDY61rMqYlIEaANMARAVY+HchBB5ITRqfbqCNt/5CKBiFQGGgHzPC4lW3zDWUuA3cC3qhrWP4/Pm8AjQIbHdQSSAtNFZKGI9PW6mACoAuwBhvmGUweLSAGvizqTSAkjE0ZEpCAwDrhPVQ96XU92qGq6qjbELcnfTETCekhVRLoAu1V1ode1BNhFqtoYt+323b5h8HAWAzQGPlDVRsARIKSvlUdKGHm+V4fxj+/ayjhghKqO97qeQPENkcwCOnlcSna1Bq72XWMZBVwqIp97W1L2qeo235+7gQm4of1wlgQkZToTH4sLp5AVKWHk+V4dJmu+C/5DgNWq+rrX9WSXiJQSkaK++/lwDTS/elpUNqnqY6oar6qVcX+PZqrqjR6XlS0iUsDXMINvKOsyIKy7VFV1J7BVRGr5DrUHQroRKCL2MwqFvToCTURGAu2AkiKSBDylqkO8rSrbWgM3Act911kAHlfVad6VlC3lgOG+bs4oYLSq5opW6FymDDDB/S5EDPCFqn7tbUkBcS8wwvcL+AbgFo/rOaOIaO02xhgT2iJlmM4YY0wIszAyxhjjOQsjY4wxnrMwMsYY4zkLI2OMMZ6zMDLGGOM5CyNjjDGe+39Df8LLVk1NbAAAAABJRU5ErkJggg==\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')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.1"
  },
  "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
}