Plotting.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "87eec03a366b4484b1184eb02573eee5",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.7, description='x_i mM', max=2.0, min=0.1), FloatSlider(value=0.1, d…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#load some modules for plotting and math\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "from scipy.integrate import odeint\n",
    "from ipywidgets import widgets\n",
    "from IPython.display import display\n",
    "import unyt as u\n",
    "plt.style.use('seaborn-poster')\n",
    "\n",
    "##Biochemical Model Function\n",
    "def BioChem_FBmodel(y,t):\n",
    "    x=y[0]\n",
    "    s=y[1]\n",
    "    V=y[2]*u.dm**3\n",
    "    ## Model parameters\n",
    "    qSmax   = 1.4   /u.hr #maximum growth rate\n",
    "    Ks      = 0.1         #saturation or Monode constant\n",
    "    qm      = 0.04  /u.hr #death rate\n",
    "    Yem     = 0.5   *u.g/u.g  #yield\n",
    "    Si      = 500         #Substrate initial\n",
    "    F0      = 0.03  *u.dm**3/u.hr #initial feed rate\n",
    "    SFR     = 0.2   *u.hr #rate of the exponential feed\n",
    "    Fmax    = 1.5   *u.dm**3/u.hr #maximum possible feed\n",
    "    \n",
    "    ## Exponential Feed Schedule\n",
    "    F = F0*numpy.exp(SFR*t)     #dm**3/hr\n",
    "    if F > Fmax:\n",
    "        F = Fmax\n",
    "    \n",
    "    qS   = qSmax*s/(s+Ks) #growth rate of your cells at given substrate concentration (2. Monod eq.) /hr\n",
    "    My   = (qS-qm)*Yem          #biomass change /hr\n",
    "\n",
    "    dy = numpy.zeros(3)\n",
    "    dy[0] = -F/V*x+My*x       \n",
    "    dy[1] = F/V *(Si-s)-qS*x   # /hr *  g/dm**3\" (dimensions \"(mass)/(length)**3\") and \"dimensionless\"\n",
    "    dy[2] = F\n",
    "    return dy\n",
    "\n",
    "def f(x_i, s_i, V_i):\n",
    "    ## Integration Time Scale\n",
    "    tspan = numpy.linspace(0, 50, 300) *u.hr\n",
    "    ## Initial Values\n",
    "    #x_i   = 0.7           #mM\n",
    "    #s_i   = 0.1           #mM\n",
    "    #V_i   = 40  *u.dm**3  #L\n",
    "    y   = [x_i,s_i,V_i]\n",
    "\n",
    "    ## Call ODE Solver\n",
    "    sol = odeint(BioChem_FBmodel,y,tspan)\n",
    "    x, s, v = sol.transpose()\n",
    "    V=v*u.dm**3  #only to get the units right, does not work with integration of function\n",
    "\n",
    "    ## Growth and Dillution Rate \n",
    "    mu=numpy.gradient(x)/(tspan[1]-tspan[0])/x                          #/hr\n",
    "    dil=numpy.gradient(V)/(tspan[1]-tspan[0])/V                         #/hr\n",
    "    X_total=numpy.array([x[i]/V[i] for i in range(x.shape[0])])         #/hr\n",
    "    deltaDilX_total=numpy.gradient(X_total)/(tspan[1]-tspan[0])/X_total\n",
    "    print(mu[1])\n",
    "\n",
    "    ## plotting\n",
    "    ##simple plot\n",
    "    #plt.plot(tspan, sol)\n",
    "    #plt.title('Volume, Substrate and Biomass')\n",
    "    #plt.xlabel('Time [hr]')\n",
    "    #plt.ylabel('Volume [liter], Concentration [mM]')\n",
    "    #plt.show()\n",
    "    ##plot of Volume, Dilution and Growth Rate\n",
    "    fig, ax1 = plt.subplots()\n",
    "    color = 'tab:red'\n",
    "    ax1.set_xlabel('Time [hr]')\n",
    "    ax1.set_ylabel('Volume, [liter]', color=color)\n",
    "    ax1.plot(tspan, V, color=color)\n",
    "    ax1.tick_params(axis='y', labelcolor=color)\n",
    "    # instantiate a second axes that shares the same x-axis\n",
    "    ax2 = ax1.twinx()  \n",
    "    color = 'tab:blue'\n",
    "    ax2.set_ylabel('Dilution and Growth Rate $\\mu$ [/hr]', color=color) \n",
    "    ax2.plot(tspan, mu, color=color, label='growth')\n",
    "    ax2.plot(tspan, dil, color=color, linestyle='--', label='dilution')\n",
    "    ax2.plot(tspan, deltaDilX_total, linestyle=':', label='difference')\n",
    "    ax2.tick_params(axis='y', labelcolor=color)\n",
    "    handles, labels = ax2.get_legend_handles_labels()\n",
    "    ax2.legend(handles, labels)\n",
    "    fig.tight_layout()  # otherwise the right y-label is slightly clipped\n",
    "    plt.show()\n",
    "    ##second plot for Biomass and Substrate\n",
    "    fig, ax1 = plt.subplots()\n",
    "    color = 'tab:red'\n",
    "    ax1.set_xlabel('Time [hr]')\n",
    "    ax1.set_ylabel('Biomass Concentration, [mM/liter]', color=color)\n",
    "    ax1.plot(tspan, x, color=color)\n",
    "    ax1.tick_params(axis='y', labelcolor=color)\n",
    "    # instantiate a second axes that shares the same x-axis\n",
    "    ax2 = ax1.twinx()  \n",
    "    color = 'tab:blue'\n",
    "    ax2.set_ylabel('Substrate Concentration [mM]', color=color) \n",
    "    ax2.plot(tspan, s, color=color)\n",
    "    ax2.tick_params(axis='y', labelcolor=color)\n",
    "    fig.tight_layout()  # otherwise the right y-label is slightly clipped\n",
    "    plt.show()\n",
    "\n",
    "interactive_plot = widgets.interact_manual(f, \n",
    "                                           x_i=widgets.FloatSlider(min=0.1, max=2, value=0.7, description='x_i mM'),\n",
    "                                           s_i=widgets.FloatSlider(min=0.01, max=0.5, value=0.1, description='s_i mM'), \n",
    "                                           V_i=widgets.FloatSlider(min=20, max=100, value=40, description='V_i L'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5b9a639700e142ca84ea62c47b02b2aa",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=58.4, description='Xyl. Tox. Lim.', max=150.0, min=20.0), FloatSlider(…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "##load some modules for plotting, math and for the interactive part (slider)\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy\n",
    "from scipy.integrate import odeint\n",
    "from ipywidgets import widgets\n",
    "from IPython.display import display\n",
    "import unyt as u\n",
    "\n",
    "#Biochemical Model Function\n",
    "def BioChem_Bmodel(x,t,Ks_d):\n",
    "    # Model Parameters Values\n",
    "    MWt_Xyl = 106.2 *u.g/u.mol    # Molecular Weight of Xylene\n",
    "    Ks_l    = 0.01                # Ks for our cells grown on Xylene (mM)\n",
    "    max_d   = 0.25  /u.hr         # Maximum Theoretical Deathrate (h-1) \n",
    "    max_g   = 0.008 /u.hr         # Maximum Theoretical Growthrate (h-1)\n",
    "    Yxx     = 2.5                 # Yield of Biomass on Xylene (--)\n",
    "    #Ks_d    = 58.4                # Ks for death due to Xylene Toxicity (mM)\n",
    "\n",
    "\n",
    "    # Growth Rate (2. Monod eq.)\n",
    "    growth = max_g * (x[1]) / (x[1] + Ks_l)\n",
    "    # Death Rate\n",
    "    death  = max_d * (x[1]) / (x[1] + Ks_d)\n",
    "\n",
    "    # Change of Substrate Concentrateion Based on Rate (3. Monod eq.)\n",
    "    # Use if you are having trouble with substrate < 0 values\n",
    "    # if x[1] < 0:\n",
    "    #     rxylene = 0\n",
    "    # else:\n",
    "    rxylene = (1/MWt_Xyl) * growth / Yxx\n",
    "\n",
    "    # Differential for Biomass and Xylene Concentration\n",
    "    biomass  = (growth - death) * x[0]   # biomass [=] Biomass (mg/L)\n",
    "    xxylene  = - rxylene * x[0]          # xxylene [=] s, Xylene Concentration (mM)\n",
    "    return [biomass, xxylene]\n",
    "\n",
    "##Integration Time Scale\n",
    "tspan = numpy.linspace(0,340,num=300)\n",
    "\n",
    "\n",
    "\n",
    "#function to update with call of widget defined below\n",
    "def f(Ks_d_move, Biomass):\n",
    "    ##Initial Conditions\n",
    "    Init_Cond = numpy.array([Biomass,1.0]) #Biomass (mg/L), Xylene Concentration (mM)\n",
    "    ##Call ODE Solver\n",
    "    solu = odeint(BioChem_Bmodel,Init_Cond,tspan,args=(Ks_d_move,))\n",
    "    x_bio, x_xyl = solu.transpose() #Xylene Concentration (mM), Biomass (mg/L)\n",
    "    ##plotting\n",
    "    fig, ax1 = plt.subplots(figsize=(14,10))\n",
    "    color = 'tab:red'\n",
    "    ax1.set_xlabel('Time [hr]')\n",
    "    ax1.set_ylabel('Cell Concentration [mg/liter]', color=color)\n",
    "    ax1.plot(tspan, x_bio, color=color)\n",
    "    #Analytical solution for unlimited exponential behavior\n",
    "    initial_growth   =  (x_bio[1]-x_bio[0])/Init_Cond[0]/(tspan[1]-tspan[0])\n",
    "    plt.ylim(top=x_bio[-1]+10)\n",
    "    ax1.plot(tspan, Init_Cond[0]*numpy.exp(initial_growth*tspan), color='red', linestyle=':', label='Exponential Growth at Initial Rate')\n",
    "    handles, labels = ax1.get_legend_handles_labels()\n",
    "    ax1.legend(handles, labels)\n",
    "    ax1.tick_params(axis='y', labelcolor=color)\n",
    "    #instantiate a second axes that shares the same x-axis\n",
    "    ax2 = ax1.twinx()  \n",
    "    color = 'tab:blue'\n",
    "    ax2.set_ylabel('Substrate Concentration [mM]', color=color)\n",
    "    ax2.plot(tspan, x_xyl, color=color)\n",
    "    ax2.tick_params(axis='y', labelcolor=color)\n",
    "\n",
    "interactive_plot = widgets.interact_manual(f, \n",
    "                                           Ks_d_move=widgets.FloatSlider(min=20, max=150, value=58.4, description='Xyl. Tox. Lim.'),\n",
    "                                           Biomass=widgets.FloatSlider(min=20, max=150, value=82.096, description='Biomass'))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import unyt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "che311",
   "language": "python",
   "name": "che311"
  },
  "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.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}