{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Enzyme Kinetics (Part 2)\n",
    "\n",
    "We have already implemented a model in which an enzyme E aids a reaction in which a substrate S is consumed to produce a product P.\n",
    "\n",
    "\\begin{eqnarray*}\n",
    "\\frac{dS}{dt} &=& -k_1\\cdot S\\cdot E + k_{-1}\\cdot SE\\\\\n",
    "\\frac{dE}{dt} &=& -k_1S\\cdot E + k_{-1}SE + k_2\\cdot SE\\\\\n",
    "\\frac{dP}{dt} &=& k_2\\cdot SE\\\\\n",
    "\\frac{dSE}{dt} &=& k_1\\cdot S\\cdot E - k_{-1}\\cdot SE - k_2\\cdot SE\\\\\n",
    "\\end{eqnarray*}\n",
    "\n",
    "where SE is the compound of S and E.\n",
    "\n",
    "Our goal today is to rewrite the model with ODEs for S and P, but not E or ES.\n",
    "\n",
    "First, we rerun the model from Tuesday. I am supplying my code so that we are all working with the equations in the same order."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.integrate\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def enzyme_model( t, y, params ):\n",
    "    S = y[0]\n",
    "    E = y[1]\n",
    "    P = y[2]\n",
    "    SE = y[3]\n",
    "\n",
    "    k1f = params[0]\n",
    "    k1b = params[1]\n",
    "    k2 = params[2]\n",
    "\n",
    "    dydt = np.zeros( (4,) )\n",
    "    dydt[0] = -k1f*S*E + k1b*SE # S\n",
    "    dydt[1] = -k1f*S*E + k1b*SE + k2*SE #E\n",
    "    dydt[2] = k2*SE #P\n",
    "    dydt[3] = k1f*S*E - k1b*SE - k2*SE #SE\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulating when Enzyme concentration is low\n",
    "Typically, the enzyme concentration is low.\n",
    "\n",
    "Use k1f = 2, k1b=2, k2 = 3, S(0) = 100, E(0) = 1, SE(0) = 0, P(0) = 0.\n",
    "\n",
    "Simulate with from t=0 to t=100 with a timestep of 0.1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k1f = 2\n",
    "k1b = 1\n",
    "k2 = 3\n",
    "params = (k1f, k1b, k2)\n",
    "\n",
    "S0 = 100\n",
    "E0 = 1\n",
    "SE0 = 0\n",
    "P0 = 0\n",
    "yinit = (S0,E0,P0,SE0)\n",
    "\n",
    "em = lambda t,y: enzyme_model(t,y,params)\n",
    "sol = scipy.integrate.solve_ivp( fun=em, t_span=(0,100), y0=yinit, \n",
    "                                 t_eval=np.arange(0,100,1) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the simulation over time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "plt.plot( sol.t, sol.y.T )\n",
    "plt.xlabel( \"Time\")\n",
    "plt.ylabel( \"Concentration\" )\n",
    "plt.legend( [\"S\",\"E\",\"P\",\"SE\"])\n",
    "plt.title( \"Enzyme-catalyzed reaction\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we zoom in on the ES complex, we see that it is nearly flat. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure()\n",
    "plt.plot( sol.t, sol.y[1,:].T )\n",
    "plt.plot( sol.t, sol.y[3,:].T )\n",
    "plt.xlabel( \"Time\")\n",
    "plt.ylabel( \"Concentration\" )\n",
    "plt.legend( [\"E\",\"SE\"])\n",
    "plt.title( \"Enzyme-catalyzed reaction\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are interested in learning how to express dP/dt in terms of S and have learned that the Michaelis-Menton equation can do that for us. If we make the assumption that it is essentially at steady-state (called the Quasi-Steady-State Assumption), then we can work out a new expression for $\\frac{dP}{dt}$\n",
    "\n",
    "The curve dP/dt vs S can be expressed as\n",
    "\n",
    "$v_{max}\\frac{S}{K_m+S}$\n",
    "\n",
    "where\n",
    "\n",
    "$v_{max} = E(0)\\cdot k_2$\n",
    "\n",
    "and \n",
    "\n",
    "$K_m = \\frac{k_{-1}+k_2}{k_1}$\n",
    "\n",
    "Let's see if that expression fits with the full 4-ODE simulation.\n",
    "\n",
    "Plot the analytic curve (from the Michaelis-Mention model for enzyme kinetics) on the same plot as the approximated dP/dt (to approximate dP/dt, compute the change in P divided by the timestep using np.diff)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Estimate the curve dP/dt from the simulation output\n",
    "P = sol.y[2,:].T\n",
    "S = sol.y[0,:-1].T\n",
    "velocity_of_P = # Write code here\n",
    "\n",
    "# Now compute it according to the Michaelis-Mention equation\n",
    "# Write code here\n",
    "analytic_dPdt = # Write code here\n",
    "plt.plot(S, analytic_dPdt)\n",
    "plt.plot( S, velocity_of_P)\n",
    "plt.xlabel( \"S\")\n",
    "plt.ylabel(\"dP/dt\")\n",
    "plt.legend( [\"M-M expression\",\"Simulation\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Model Enzyme Kinetics with the Michaelis-Menton Equation\n",
    "\n",
    "Implement a model that uses the M-M equation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the same initial conditions as above, and computing Km and vmax according to\n",
    "\n",
    "$v_{max} = E(0)\\cdot k_2$\n",
    "\n",
    "and \n",
    "\n",
    "$K_m = \\frac{k_{-1}+k_2}{k_1}$\n",
    "\n",
    "simulate the M-M model for the same time steps as the full model."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the new model's simulation results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now compare P from one model to the other by plotting them both"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Awesome, the full model and the M-M model are nearly indistinguishable. When we are modeling a process using M-M kinetics, we need two parameters, $v_{max}$ and $K_m$. Above, we saw that the two parameters depend on a larger set of parameters and the initial concentration of E. Typically, we will be assigning values directly to $v_{max}$ and $K_m$, skipping the step of determining separate parameters for forward and backward reaction rates. We should develop and intuition for how $v_{max}$ and $K_m$ affect dP/dt and P(t).\n",
    "\n",
    "$v_{max}$ is simply the maximum value of dP/dt.\n",
    "\n",
    "$K_m$ is often referred to as the activation threshold because it identifies the concentration of S that causes dP/dt to be at half it maximum. \n",
    "\n",
    "Make an array of S values and compute the M-M equationn for a range of values of Km. The trend is that curves with larger Km are less steep."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copy-paste from notebook with code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the simulation with a small Km and compare it to a simulation with a large Km"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Copy-paste from notebook with code or write it on your own after class"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Compare to S -> P\n",
    "\n",
    "This is bonus material that shows that enzyme-aided production is fundamentally different from simple production. The rate of production can become saturated (reach an upper limit) with enyme-aided production, but not with simple production. We see this by plotting the rate of production of P in the simple production model. The rate is linear. This is in contrast to the similar plot from enyme-aided production, in which the rate reaches a threshold and flattens out."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def production_model( t, y, params ):\n",
    "    S = y[0]\n",
    "    P = y[1]\n",
    "\n",
    "    k1 = params[0]\n",
    "\n",
    "    dydt = np.zeros( (2,) )\n",
    "    dydt[0] = -k1*S # S\n",
    "    dydt[1] = k1*S #P\n",
    "    return dydt\n",
    "\n",
    "k1 = 0.1\n",
    "params = (k1,)\n",
    "\n",
    "S0 = 100\n",
    "P0 = 0\n",
    "yinit = (S0,P0)\n",
    "\n",
    "em = lambda t,y: production_model(t,y,params)\n",
    "sol = scipy.integrate.solve_ivp( fun=em, t_span=(0,100), y0=yinit, \n",
    "                                 t_eval=np.arange(0,100,1) )\n",
    "\n",
    "plt.plot( sol.t, sol.y.T )\n",
    "plt.xlabel( \"Time\")\n",
    "plt.ylabel( \"Concentration\" )\n",
    "plt.legend( [\"S\",\"P\"])\n",
    "plt.title( \"Production reaction\");\n",
    "\n",
    "# Estimate the curve dP/dt from the simulation output\n",
    "P = sol.y[1,:].T\n",
    "S = sol.y[0,:-1].T\n",
    "velocity_of_P = np.divide(np.diff(P), (sol.t[1]-sol.t[0]))\n",
    "RHS = k1*S\n",
    "\n",
    "# Now compute it according to the Michaelis-Mention equation\n",
    "plt.figure()\n",
    "plt.plot(S, velocity_of_P)\n",
    "plt.plot(S, RHS)\n",
    "plt.xlabel( \"S\")\n",
    "plt.ylabel(\"dP/dt\")\n",
    "plt.legend( [\"Simulation\",\"RHS\"]);\n",
    "plt.title( \"Rate of Production Reaction\")"
   ]
  }
 ],
 "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}