{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "bc848aef-c031-4db3-bf24-61bd0e4cf0fe",
   "metadata": {},
   "source": [
    "# Michaelis-Menten"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4427cead-5930-449c-a2b9-64a79a8ec6a1",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "from pykappa import System"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "71ff78ff-d429-4c0f-bd3f-c7d9c2ec0274",
   "metadata": {},
   "source": [
    "In this example, which is adapted from the [Kappa language manual](https://kappalanguage.org/static/manual.pdf), we consider the standard Michaelis-Menten system and an approximation. The first model is precise in that an enzyme `E` binds a substrate `S` and subsequently phosphorylates it (site `x{u}` becomes `x{p}`). This model explicitly specifies a binding requirement for phosphorylation. The second model, in contrast, collapses the binding and phosphorylation step into a single rule. The dynamics are instead specified by the rate term, which is based on the quasi steady-state approximation for Michaelis-Menten kinetics, where the rate of phosphorylation is\n",
    "\n",
    "$$ \\frac{k_2 E_t [S]}{K_m + [S]}, $$\n",
    "\n",
    "where $k_2$ is a catalytic rate constant, $[S]$ is the concentration of free substrate, and $E_t$ is the total enzyme concentration."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0eeea96e-0c48-43cd-8209-7fc05360b7cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "kastr = \"\"\"\n",
    "    %var: 'k1' 0.001\n",
    "    %var: 'k_1' 0.1\n",
    "    %var: 'k2' 1.0\n",
    "    %var: 'Km' ('k_1'+'k2')/'k1'\n",
    "\n",
    "    // Precise model\n",
    "    \n",
    "    E(s[.]), S(e[.],x{u}) <-> E(s[1]), S(e[1],x{u}) @ 'k1', 'k_1'\n",
    "    E(s[1]), S(e[1],x{u}) -> E(s[.]), S(e[.],x{p})  @ 'k2'\n",
    "\n",
    "    %obs: 'P' |S(e[.],x{p})|\n",
    "    %obs: 'S' |S(e[.], x{u})|\n",
    "    %obs: 'ES' |E(s[1]), S(e[1],x{u})|\n",
    "\n",
    "    // Approximate model\n",
    "    \n",
    "    _E(s[.]), _S(e[.],x{u}) -> _E(s[.]), _S(e[.],x{p}) @ 'k2'/('Km'+'~S')\n",
    "\n",
    "    %obs: '~S' |_S(e[.],x{u})|\n",
    "    %obs: '~P' |_S(e[.],x{p})|\n",
    "\"\"\""
   ]
  },
  {
   "cell_type": "markdown",
   "id": "560af8ac-d5c1-4e3f-9325-5c161234221e",
   "metadata": {},
   "source": [
    "Note that $E_t [S]$ has been omitted from the Kappa model since PyKappa already applies the rule at a rate exactly proportional to the number of embeddings $E_t [S]$. (All enzymes are free in the approximate model.)\n",
    "\n",
    "We now initialize this model with an excess of substrate, a requirement for the quasi steady-state approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4c37e8f8-3c8d-4c3b-91a7-45ce1a549cbb",
   "metadata": {},
   "outputs": [],
   "source": [
    "qe_system = System.from_ka(kastr, seed=42)\n",
    "\n",
    "qe_system.mixture.add(\"E(s[.])\", 100)\n",
    "qe_system.mixture.add(\"S(e[.], x{u})\", 100000)\n",
    "qe_system.mixture.add(\"_E(s[.])\", 100)\n",
    "qe_system.mixture.add(\"_S(e[.], x{u})\", 100000)\n",
    "\n",
    "while qe_system.reactivity:\n",
    "    qe_system.update()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "fd78ef8d-f88e-4320-911f-d73abb9884e8",
   "metadata": {},
   "source": [
    "Plotting these two models against each other demonstrates the reliability of the quasi steady-state approximation under appropriate conditions. Note that we add a separate axis for the complex `ES`; the scale differs from the other quantities by several orders of magnitude."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1125637c-57c0-491a-a64d-1cd37d3aa0b5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fig = qe_system.monitor.plot(combined=True, observables=[\"P\", \"S\", \"~P\"]);\n",
    "\n",
    "ax = fig.axes[0]\n",
    "ax_twin = ax.twinx()\n",
    "ax_twin.margins(0, 0)\n",
    "ax_twin.plot(qe_system.monitor.history[\"time\"], qe_system.monitor.history[\"ES\"], \n",
    "             color=\"#d62728\", label=\"ES\", linewidth=0.5)\n",
    "lines1, labels1 = ax.get_legend_handles_labels()\n",
    "lines2, labels2 = ax_twin.get_legend_handles_labels()\n",
    "ax.legend(lines1 + lines2, labels1 + labels2, loc=\"center left\");"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "bee6b6d8-42ab-48bb-b22a-8e98c303743f",
   "metadata": {},
   "source": [
    "Note that the green line (`~P`, the number of phosphorylated agents under a less detailed ruleset) closely tracks the blue line (`P`, the higher-fidelity model). This model breaks down when the quasi steady-state dynamics are not satisfied, which can occur under an overabundance of enzymes compared to substrates: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e2a5f58e-e1fb-4b80-bcb8-67858f82df10",
   "metadata": {},
   "outputs": [],
   "source": [
    "nqe_system = System.from_ka(kastr, seed=42)\n",
    "\n",
    "# disrupt the quasi-equilibrium\n",
    "nqe_system[\"k2\"] = \"0.1\"\n",
    "\n",
    "nqe_system.mixture.add(\"E(s[.])\", 10000)\n",
    "nqe_system.mixture.add(\"S(e[.], x{u})\", 1000)\n",
    "nqe_system.mixture.add(\"_E(s[.])\", 10000)\n",
    "nqe_system.mixture.add(\"_S(e[.], x{u})\", 1000)\n",
    "\n",
    "while nqe_system.time < 25:\n",
    "    nqe_system.update()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43eaf1ae-1fe4-451f-b71c-e006364fc9dd",
   "metadata": {},
   "outputs": [],
   "source": [
    "nqe_system.monitor.plot(combined=True, observables=[\"P\", \"S\", \"~P\", \"ES\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "0b773fdf-2603-4976-9659-31d3b1217c9d",
   "metadata": {},
   "source": [
    "Note that `~P` no longer tracks the more correct estimate `P`."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "standard",
   "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.14.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}