{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from scipy.special import binom"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Independent binding sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 22196 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 28900 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 65424 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "L = 10**(np.linspace(-5, -1, 100))\n",
    "K = 1000\n",
    "n = 5\n",
    "\n",
    "Q = (1 + K * L)**n\n",
    "\n",
    "P = np.zeros((n + 1, L.size))\n",
    "P[0] = 1 / Q\n",
    "for i in range(1, n+1):\n",
    "    P[i] = binom(n, i) * K**i * L**i / Q\n",
    "\n",
    "nu = np.zeros_like(L)\n",
    "for i in range(1, n+1):\n",
    "    nu += i * binom(n, i) * K**i * L**i / Q\n",
    "\n",
    "\n",
    "plt.title(\"Bound ligands\")\n",
    "plt.xlabel(\"[L]\")\n",
    "plt.xscale(\"log\")\n",
    "plt.ylabel(\"$\\\\nu$\")\n",
    "plt.plot(L, nu)\n",
    "plt.show()\n",
    "\n",
    "plt.title(\"Scatchard Plot\")\n",
    "plt.xlabel(\"$\\\\nu$\")\n",
    "plt.ylabel(\"$\\\\nu$/[L]\")\n",
    "plt.plot(nu, nu / L)\n",
    "plt.show()\n",
    "\n",
    "\n",
    "plt.title(\"Probability of single states\")\n",
    "plt.xlabel(\"[L]\")\n",
    "plt.xscale('log')\n",
    "for i in range(0, n+1):\n",
    "    plt.plot(L, P[i], label=f\"{i}\")\n",
    "plt.legend()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two kinds of independent binding sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 36904 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 38732 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 74356 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "L = 10**(np.linspace(-5, -1, 100))\n",
    "K_A = 1000\n",
    "n_A = 5\n",
    "K_B = 10000\n",
    "n_B = 2\n",
    "\n",
    "\n",
    "Q_A = (1 + K_A * L)**n_A\n",
    "\n",
    "P_A = np.zeros((n_A + 1, L.size))\n",
    "P_A[0] = 1 / Q_A\n",
    "for i in range(1, n_A+1):\n",
    "    P_A[i] = binom(n_A, i) * K_A**i * L**i / Q_A\n",
    "\n",
    "nu_A = np.zeros_like(L)\n",
    "for i in range(1, n_A+1):\n",
    "    nu_A += i * binom(n_A, i) * K_A**i * L**i / Q_A\n",
    "\n",
    "Q_B = (1 + K_B * L)**n_B\n",
    "\n",
    "P_B = np.zeros((n_B + 1, L.size))\n",
    "P_B[0] = 1 / Q_B\n",
    "for i in range(1, n_B+1):\n",
    "    P_B[i] = binom(n_B, i) * K_B**i * L**i / Q_B\n",
    "\n",
    "nu_B = np.zeros_like(L)\n",
    "for i in range(1, n_B+1):\n",
    "    nu_B += i * binom(n_B, i) * K_B**i * L**i / Q_B\n",
    "\n",
    "\n",
    "nu = nu_A + nu_B\n",
    "\n",
    "P = np.zeros((n_A + n_B + 1, L.size))\n",
    "for i in range(0, n_A + n_B + 1):\n",
    "    for j in range(0,i+1):\n",
    "        if j < n_A + 1 and i-j < n_B + 1:\n",
    "            P[i] += P_A[j] * P_B[i-j] \n",
    "           \n",
    "\n",
    "plt.title(\"Bound ligands\")\n",
    "plt.xlabel(\"[L]\")\n",
    "plt.xscale(\"log\")\n",
    "plt.ylabel(\"$\\\\nu$\")\n",
    "plt.plot(L, nu, label=\"$\\\\nu$\")\n",
    "plt.plot(L, nu_A, label=\"$\\\\nu_A$\")\n",
    "plt.plot(L, nu_B, label=\"$\\\\nu_B$\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.title(\"Scatchard Plot\")\n",
    "plt.xlabel(\"$\\\\nu$\")\n",
    "plt.ylabel(\"$\\\\nu$/[L]\")\n",
    "plt.plot(nu, nu / L, label=\"$\\\\nu / L$\")\n",
    "plt.plot(nu_A, nu_A / L, label=\"$\\\\nu_A / L$\")\n",
    "plt.plot(nu_B, nu_B / L, label=\"$\\\\nu_B / L$\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "\n",
    "plt.title(\"Probability of single states\")\n",
    "plt.xlabel(\"[L]\")\n",
    "plt.xscale('log')\n",
    "for i in range(0, n_A + n_B + 1):\n",
    "    plt.plot(L, P[i], label=f\"{i}\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Infinitely cooperative binding sites"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 23252 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 29796 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "L = 10**(np.linspace(-5, -1, 100))\n",
    "K = 1000\n",
    "n = 5\n",
    "\n",
    "Q = np.ones_like(L) + K**n * L**n\n",
    "P = K**n * L**n / Q\n",
    "\n",
    "\n",
    "plt.title(\"Probability for occupied binding sites\")\n",
    "plt.xlabel(\"[L]\")\n",
    "plt.xscale('log')\n",
    "plt.plot(L, P)\n",
    "plt.show()\n",
    "\n",
    "plt.title(\"Hill-Plot\")\n",
    "plt.xlabel(\"[L]\")\n",
    "plt.xscale('log')\n",
    "plt.yscale('log')\n",
    "plt.ylabel(\"$\\\\log\\\\left(\\\\frac{ P }{1-P}\\\\right) = n (\\\\log L + \\\\log K)$\")\n",
    "plt.plot(L, P / (1-P))\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hemoglobin"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "[encoded data, 26180 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 47400 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 24972 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 49044 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "[encoded data, 28252 characters]",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "L = 10**(np.linspace(-12, -1, 100))\n",
    "K0 = 10e-4\n",
    "K1 = 10e5\n",
    "K2 = 10e7\n",
    "n = 4\n",
    "\n",
    "Q1 = (1 + L * K1)**n\n",
    "Q2 = (1 + L * K2)**n\n",
    "\n",
    "Qtot = Q1 + K0 * Q2\n",
    "\n",
    "P1 = np.zeros((n + 1, L.size))\n",
    "P1[0] = 1 / Qtot\n",
    "for i in range(1, n+1):\n",
    "    P1[i] = binom(n, i) * K1**i * L**i / Qtot\n",
    "P1_tot = np.sum(P1, axis=0)\n",
    "\n",
    "P2 = np.zeros((n + 1, L.size))\n",
    "P2[0] = K0 / Qtot\n",
    "for i in range(1, n+1):\n",
    "    P2[i] = K0 * binom(n, i) * K2**i * L**i / Qtot\n",
    "P2_tot = np.sum(P2, axis=0)\n",
    "\n",
    "nu1 = np.zeros_like(L)\n",
    "for i in range(1, n+1):\n",
    "    nu1 += i * binom(n, i) * K1**i * L**i / Qtot\n",
    "\n",
    "nu2 = np.zeros_like(L)\n",
    "for i in range(1, n+1):\n",
    "    nu2 += i * binom(n, i) * K0 * K2**i * L**i / Qtot\n",
    "\n",
    "nu = nu1 + nu2\n",
    "\n",
    "P = nu / n\n",
    "\n",
    "plt.title(\"Bound oxygen\")\n",
    "plt.xlabel(\"[L]\")\n",
    "plt.xscale(\"log\")\n",
    "plt.ylabel(\"$\\\\nu$\")\n",
    "plt.plot(L, nu, label=\"$\\\\nu_{tot}$\")\n",
    "plt.plot(L, nu1, label=\"$\\\\nu_1$\")\n",
    "plt.plot(L, nu2, label=\"$\\\\nu_2$\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.title(\"Scatchard Plot\")\n",
    "plt.xlabel(\"$\\\\nu$\")\n",
    "plt.ylabel(\"$\\\\nu$/[L]\")\n",
    "plt.plot(nu, nu / L, label=\"$\\\\nu_{tot} / L$\")\n",
    "plt.plot(nu1, nu1 / L, label=\"$\\\\nu_1 / L$\")\n",
    "plt.plot(nu2, nu2 / L, label=\"$\\\\nu_2 / L$\")\n",
    "plt.legend()\n",
    "plt.show() \n",
    "\n",
    "plt.title(\"Hill Plot\")\n",
    "plt.xlabel(\"[L]\")\n",
    "plt.xscale('log')\n",
    "plt.ylabel(\"$\\\\frac{P}{1-P} = \\\\frac{\\\\nu}{n-\\\\nu}$\")\n",
    "plt.yscale('log')\n",
    "plt.plot(L, P / (1-P))\n",
    "plt.show()\n",
    "\n",
    "\n",
    "plt.title(\"Probability of single states\")\n",
    "plt.xlabel(\"[L]\")\n",
    "plt.xscale('log')\n",
    "for i in range(0, n+1):\n",
    "    plt.plot(L, P1[i], label=f\"Deoxy {i}\")\n",
    "for i in range(0, n+1):\n",
    "    plt.plot(L, P2[i], label=f\"Oxy {i}\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "plt.title(\"Probability of Deoxy vs Oxy\")\n",
    "plt.xlabel(\"[L]\")\n",
    "plt.xscale('log')\n",
    "plt.plot(L, P1_tot, label=f\"Deoxy\")\n",
    "plt.plot(L, P2_tot, label=f\"Oxy\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "venv",
   "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.10"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
