test_basic_functions.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "test_basic_functions.ipynb",
      "version": "0.3.2",
      "provenance": [],
      "toc_visible": true,
      "include_colab_link": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/github/fbeilstein/machine_learning/blob/master/test_basic_functions.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "metadata": {
        "id": "vipfqylf8pXk",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "---\n",
        "\n",
        "To test plotting I've added the code that I've written as practice during my study at the \"Biology And Medicine Institute\" (Science Educational Center Of Taras Shevchenko National University Of Kyiv).\n",
        "\n",
        "---"
      ]
    },
    {
      "metadata": {
        "id": "YHFw3qmJJJMu",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **1.** The ligand-binding equation\n",
        "\n",
        "For more details please see Wikipedia:\n",
        "\n",
        "https://en.wikipedia.org/wiki/Michaelis%E2%80%93Menten_kinetics\n",
        "\n",
        "https://en.wikipedia.org/wiki/Hill_equation_(biochemistry)\n",
        "\n"
      ]
    },
    {
      "metadata": {
        "id": "bKt8Omb5LE-g",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "# import packages needed for calculation and plotting\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "# In colabs this lines cause ERROR\n",
        "# Nevertheless, LaTeX works (but badly)\n",
        "#plt.rcParams['text.usetex'] = True\n",
        "#plt.rcParams['text.latex.unicode'] = True\n",
        "\n",
        "# start figure\n",
        "plt.figure(1, figsize=(8, 5))\n",
        "# position axes origin and (xmax;ymax)\n",
        "ax = plt.axes([0.11, 0.11, 0.8, 0.75])\n",
        "# create arrays of arguments and parameters\n",
        "L = np.arange(0.0, 200.0, 0.1)\n",
        "Kd = np.arange(10.0, 91.0, 20.0)\n",
        "Rtot = 20 * 10**-3\n",
        "# generate array of values for this arguments\n",
        "Lr = [Rtot * L / (_Kd + L) for _Kd in Kd]\n",
        "# perform plotting and generate legends\n",
        "plots = [plt.plot(L, _Lr)[0] for _Lr in Lr]\n",
        "legends = [r\"$K_D = \" + str(_Kd)\\\n",
        "                      + r\"~\\mu M$\" for _Kd in Kd]\n",
        "plt.legend(plots, legends, loc=4)\n",
        "\n",
        "# set additional parameters such as labels\n",
        "plt.xlabel(r'$[L], (\\mu M)$', fontsize=16)\n",
        "plt.ylabel(r'$[LR], (nM)$', fontsize=16)\n",
        "# \\displaystyle and \\cfrac cause error\n",
        "# LaTeX is patially supported\n",
        "plt.title(r'Ligand-Receptor binding '+\\\n",
        "   r'$[LR] = \\frac{R_{TOT}[L]}{K_D + [L]}$',\\\n",
        "   fontsize=16, color='black')\n",
        "plt.grid(True)\n",
        "# save for later usage\n",
        "plt.draw()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "TrtY3R-RMYEk",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **2.** Lorentz attractor\n",
        "\n",
        "For more details please see Wikipedia:\n",
        "\n",
        "https://en.wikipedia.org/wiki/Lorenz_system#Lorenz_attractor"
      ]
    },
    {
      "metadata": {
        "id": "bsJNFketNT6O",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "from mpl_toolkits.mplot3d import Axes3D\n",
        "\n",
        "\n",
        "# define function calculating derivatives\n",
        "def lorenz(variables, t, sigma=10.0):\n",
        "    sigma = 10.0\n",
        "    rho = 28.0\n",
        "    beta = 8.0/3.0\n",
        "    [x, y, z] = variables\n",
        "    # calculate derivatives\n",
        "    xd = sigma * (y-x)\n",
        "    yd = (rho-z)*x - y\n",
        "    zd = x*y - beta*z\n",
        "    # return the states\n",
        "    return [xd, yd, zd]\n",
        "\n",
        "initial_position = [2.0, 3.0, 4.0]\n",
        "t = np.arange(0.0, 30.0, 0.01)\n",
        "solution = odeint(lorenz, initial_position, t)\n",
        "result = np.transpose(solution)\n",
        "\n",
        "# plot the result\n",
        "fig = plt.figure()\n",
        "ax = fig.gca(projection='3d')\n",
        "ax.plot(result[0], result[1], result[2])\n",
        "ax.set_xlabel('x')\n",
        "ax.set_ylabel('y')\n",
        "ax.set_zlabel('z')\n",
        "# save results\n",
        "plt.draw()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "sQNzwnXnQ2ai",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **3.** Myxococcus xantus C-signal pathway\n",
        "\n",
        "For more details please see article: \n",
        "\n",
        "Mogilner, Alex, Roy Wollman, and Wallace F. Marshall. \"Quantitative modeling in cell biology: what is it good for?.\" *Developmental cell* 11, no. 3 (2006): 279-287."
      ]
    },
    {
      "metadata": {
        "id": "yzzu1QspRCNC",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "\n",
        "# define function calculating derivatives\n",
        "def func_Frz(variables, t):\n",
        "    k = [1, 4, 4]\n",
        "    k_bar = [0.08, 2, 2]\n",
        "    K = [0.005, 0.005, 0.005]\n",
        "    K_bar = [0.01, 0.005, 0.005]\n",
        "    [Frz, FrzCD, FrzE] = variables\n",
        "    # calculate derivatives\n",
        "    dFrz = k[0] * (1-FrzE) * (1-Frz) / (1-Frz+K[0])\\\n",
        "         - k_bar[0] * Frz / (Frz+K_bar[0])\n",
        "    dFrzCD = k_bar[1] * (1-FrzCD) / (1-FrzCD+K_bar[1])\\\n",
        "           - k[1] * (1-Frz) * FrzCD / (FrzCD+K[1])\n",
        "    dFrzE = k_bar[2] * (1-FrzE) / (1-FrzE+K_bar[2])\\\n",
        "          - k[2] * (1-FrzCD) * FrzE / (FrzE+K[2])\n",
        "    # return the states\n",
        "    return [dFrz, dFrzCD, dFrzE]\n",
        "\n",
        "initial_position = [1, 1, 1]\n",
        "t = np.arange(0.0, 30.0, 0.01)\n",
        "solution = odeint(func_Frz, initial_position, t)\n",
        "result = np.transpose(solution)\n",
        "\n",
        "plots = [plt.plot(t, _result)[0] for _result in result]\n",
        "legends = [r'[Frz]/[Frz]tot', r'[FrzCD]/[FrzCD]tot',\\\n",
        "           r'[FrzE]/[FrzE]tot']\n",
        "plt.legend(plots, legends, loc=1)\n",
        "plt.xlabel(r'time, $t$, min', fontsize=16)\n",
        "plt.ylabel(r'relative concentration, [X]/[X]tot',\\\n",
        "           fontsize=16)\n",
        "# save results\n",
        "plt.draw()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "JutG6HVuTBRk",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **4.** Myxococcus xantus C-signal pathway and external influence\n",
        "\n",
        "For more details please see article:\n",
        "\n",
        "Igoshin, Oleg A., Albert Goldbeter, Dale Kaiser, and George Oster. \"A biochemical oscillator explains several aspects of Myxococcus xanthus behavior during development.\" *Proceedings of the National Academy of Sciences* 101, no. 44 (2004): 15760-15765."
      ]
    },
    {
      "metadata": {
        "id": "ON_j84_fTXxi",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "\n",
        "# define function calculating derivatives\n",
        "def func_Frz(variables, t):\n",
        "    k = [1, 4, 4]\n",
        "    k_bar = [0.08, 2, 2]\n",
        "    if t > 30:\n",
        "        k_bar = [0.58, 2, 2]\n",
        "    K = [0.005, 0.005, 0.005]\n",
        "    K_bar = [0.01, 0.005, 0.005]\n",
        "    [Frz, FrzCD, FrzE] = variables\n",
        "    # calculate derivatives\n",
        "    dFrz = k[0] * (1-FrzE) * (1-Frz) / (1-Frz+K[0])\\\n",
        "         - k_bar[0] * Frz / (Frz+K_bar[0])\n",
        "    dFrzCD = k_bar[1] * (1-FrzCD) / (1-FrzCD+K_bar[1])\\\n",
        "           - k[1] * (1-Frz) * FrzCD / (FrzCD+K[1])\n",
        "    dFrzE = k_bar[2] * (1-FrzE) / (1-FrzE+K_bar[2])\\\n",
        "          - k[2] * (1-FrzCD) * FrzE / (FrzE+K[2])\n",
        "    # return the states\n",
        "    return [dFrz, dFrzCD, dFrzE]\n",
        "\n",
        "initial_position = [1, 1, 1]\n",
        "t = np.arange(0.0, 60.0, 0.01)\n",
        "solution = odeint(func_Frz, initial_position, t)\n",
        "result = np.transpose(solution)\n",
        "\n",
        "plots = [plt.plot(t, _result)[0] for _result in result]\n",
        "legends = [r'[Frz]/[Frz]tot', r'[FrzCD]/[FrzCD]tot',\\\n",
        "           r'[FrzE]/[FrzE]tot']\n",
        "plt.legend(plots, legends, loc=1)\n",
        "plt.xlabel(r'time, $t$, min', fontsize=16)\n",
        "plt.ylabel(r'relative concentration, [X]/[X]tot',\\\n",
        "           fontsize=16)\n",
        "plt.plot([30, 30], [0, 1], 'k:', lw=2)\n",
        "plt.text(28, 0.2, r'SIGNAL', fontsize=20, rotation=90)\n",
        "# save results\n",
        "plt.draw()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "q73bCSC9T_6k",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **5.** Saccharomyces cerevisiae glycolytic oscillations\n",
        "\n",
        "For more details please see article:\n",
        "\n",
        "Bier, Martin, Barbara M. Bakker, and Hans V. Westerhoff. \"How yeast cells synchronize their glycolytic oscillations: a perturbation analytic treatment.\"* Biophysical Journal* 78, no. 3 (2000): 1087-1093."
      ]
    },
    {
      "metadata": {
        "id": "iuNBFhi5UUdN",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "\n",
        "def func_G_ATP(variables, t, Km):\n",
        "    Vin = 0.36\n",
        "    k1 = 0.02\n",
        "    kp = 6\n",
        "    [G, ATP] = variables\n",
        "    # calculate derivatives\n",
        "    dG = Vin - k1 * G * ATP\n",
        "    dATP = 2 * k1 * G *ATP - kp * ATP / (ATP + Km)\n",
        "    # return the states\n",
        "    return [dG, dATP]\n",
        "\n",
        "initial_position = [1, 1]\n",
        "t = np.arange(0.0, 600.0, 1)\n",
        "K_low = 12\n",
        "K_high = 20\n",
        "low_const = lambda v, time: func_G_ATP(v, time, Km=K_low)\n",
        "high_const = lambda v, time: func_G_ATP(v, time, Km=K_high)\n",
        "\n",
        "solution = odeint(low_const, initial_position, t)\n",
        "result = np.transpose(solution)\n",
        "plots = [plt.plot(t, _result)[0] for _result in result]\n",
        "legends = [r'[G]', r'[ATP]']\n",
        "plt.legend(plots, legends, loc=1)\n",
        "plt.xlabel(r'time, $t$, min', fontsize=16)\n",
        "plt.ylabel(r'relative concentration, [X]', fontsize=16)\n",
        "plt.text(300, 21, r\"$K_m = $\" + str(K_low), fontsize=20)\n",
        "# save results\n",
        "plt.draw()\n",
        "plt.savefig('day_7_1.eps')\n",
        "\n",
        "plt.clf()\n",
        "solution = odeint(high_const, initial_position, t)\n",
        "result = np.transpose(solution)\n",
        "plots = [plt.plot(t, _result)[0] for _result in result]\n",
        "legends = [r'[G]', r'[ATP]']\n",
        "plt.legend(plots, legends, loc=1)\n",
        "plt.xlabel(r'time, $t$, min', fontsize=16)\n",
        "plt.ylabel(r'relative concentration, [X]', fontsize=16)\n",
        "plt.text(300, 12, r\"$K_m = $\" + str(K_high), fontsize=20)\n",
        "# save results\n",
        "plt.draw()\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "uqU9QTnOUzIX",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **6.** Saccharomyces cerevisiae glycolytic oscillations, phase portrait\n",
        "\n",
        "For more details please see article:\n",
        "\n",
        "Bier, Martin, Barbara M. Bakker, and Hans V. Westerhoff. \"How yeast cells synchronize their glycolytic oscillations: a perturbation analytic treatment.\"* Biophysical Journal* 78, no. 3 (2000): 1087-1093."
      ]
    },
    {
      "metadata": {
        "id": "s6-LRyeMU5nT",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "\n",
        "# define function calculating derivatives\n",
        "def func_G_ATP(variables, t, Km, Vin=0.36, k1=0.02, kp=6):\n",
        "    [G, ATP] = variables\n",
        "    dG = Vin - k1 * G * ATP\n",
        "    dATP = 2 * k1 * G *ATP - kp * ATP / (ATP + Km)\n",
        "    return [dG, dATP]\n",
        "\n",
        "# define function for painting arrows\n",
        "def calc_arrows(data_array, threshold):\n",
        "    solution_diffs = np.hstack((data_array[:-1, :],\\\n",
        "         np.transpose(np.diff(np.transpose(data_array)))))\n",
        "    arrows = []\n",
        "    for [x, y, dx, dy] in solution_diffs:\n",
        "        use_this = True\n",
        "        for [ax, ay, dax, day] in arrows:\n",
        "            if (ax - x)**2 + (ay - y)**2 < threshold:\n",
        "                use_this = False\n",
        "                break\n",
        "        if use_this:\n",
        "            norm = np.sqrt(dx**2 + dy**2)\n",
        "            arrows.append([x, y, dx/norm, dy/norm])\n",
        "    return np.transpose(arrows)\n",
        "\n",
        "# solving and plotting function\n",
        "def do_plotting(_Km, _initial, _time, _scale,\n",
        "                _filename, _tx, _ty):\n",
        "    l_func = lambda v, time: func_G_ATP(v, time, Km=_Km)\n",
        "    solution = odeint(l_func, _initial, _time)\n",
        "    plt.plot(solution[:, 0], solution[:, 1])\n",
        "    arrows = calc_arrows(solution, .5)\n",
        "    plt.quiver(arrows[0],arrows[1],arrows[2],arrows[3],\\\n",
        "               scale_units='xy',angles='xy',scale=_scale)\n",
        "    plt.xlabel(r'glucose, [G]', fontsize=16)\n",
        "    plt.ylabel(r'ATP, [ATP]', fontsize=16)\n",
        "    plt.text(_tx, _ty, r\"$K_m = \" + str(_Km) + \"$\",\\\n",
        "             fontsize=20)\n",
        "    plt.draw()\n",
        "    plt.savefig(_filename)\n",
        "\n",
        "initial_position = [10, 5]\n",
        "t = np.arange(0.0, 600.0, 0.1)\n",
        "K_low = 12\n",
        "K_high = 20\n",
        "\n",
        "#do_plotting(K_low, initial_position, t, 1,\\\n",
        "#            'day_8_1.eps', 20, 14)\n",
        "#plt.clf()\n",
        "#do_plotting(K_high, initial_position, t, 5,\\\n",
        "#            'day_8_2.eps', 6, 7)\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "CYtURjmDVSNN",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "do_plotting(K_low, initial_position, t, 1,\\\n",
        "            'day_8_1.eps', 20, 14)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "Oaq3TvMmVVtL",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "do_plotting(K_high, initial_position, t, 5,\\\n",
        "            'day_8_2.eps', 6, 7)"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "kf0ly7TQVsL1",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **7.** Saccharomyces cerevisiae glycolytic oscillations, bifurcation\n",
        "\n",
        "For more details please see article:\n",
        "\n",
        "Bier, Martin, Barbara M. Bakker, and Hans V. Westerhoff. \"How yeast cells synchronize their glycolytic oscillations: a perturbation analytic treatment.\"* Biophysical Journal* 78, no. 3 (2000): 1087-1093."
      ]
    },
    {
      "metadata": {
        "id": "_VaOb1MvVzBB",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "\n",
        "# define function calculating derivatives\n",
        "def func_G_ATP(variables, t, Km, Vin=0.36, k1=0.02, kp=6):\n",
        "    [G, ATP] = variables\n",
        "    # calculate derivatives\n",
        "    dG = Vin - k1 * G * ATP\n",
        "    dATP = 2 * k1 * G *ATP - kp * ATP / (ATP + Km)\n",
        "    # return the states\n",
        "    return [dG, dATP]\n",
        "\n",
        "# fixed point\n",
        "def fixed_point(_Km, Vin=0.36, k1=0.02, kp=6):\n",
        "    return (kp - 2 * Vin) / (2 * k1 * _Km)\n",
        "\n",
        "Km = np.arange(10, 25, .1)\n",
        "t = np.arange(0.0, 600.0, 0.1)\n",
        "initial_position = [10, 5]\n",
        "\n",
        "KmMAX = []\n",
        "KmMIN = []\n",
        "for _Km in Km:\n",
        "    l_func = lambda v, time: func_G_ATP(v, time, Km=_Km)\n",
        "    solution = odeint(l_func, initial_position, t)\n",
        "    KmMAX.append(np.amax(solution[5000:, 0]))\n",
        "    KmMIN.append(np.amin(solution[5000:, 0]))\n",
        "K_fixed = fixed_point(Km)\n",
        "\n",
        "plt.plot(Km, KmMAX)\n",
        "plt.plot(Km, KmMIN)\n",
        "plt.fill_between(Km, KmMIN, KmMAX, facecolor='yellow')\n",
        "plt.plot(Km, K_fixed)\n",
        "plt.ylabel(r'glucose, [G]', fontsize=16)\n",
        "plt.xlabel(r'Km', fontsize=16)\n",
        "plt.text(13, 18, r\"$[G]_{max}$\", fontsize=20)\n",
        "plt.text(14, 3, r\"$[G]_{min}$\", fontsize=20)\n",
        "plt.text(12, 12, r\"$[G]_{fixed}$\", fontsize=20)\n",
        "plt.text(11, 15, r\"Oscillations\", fontsize=20, rotation=90)\n",
        "plt.text(20, 8, r\"Fixed point\", fontsize=20)\n",
        "plt.draw()"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "dFaW5475WUkn",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **8.** Bio-switch with co-acivation\n",
        "\n",
        "For more details please see article:\n",
        "\n",
        "Tyson, John J., Katherine C. Chen, and Bela Novak. \"Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell.\" *Current opinion in cell biology* 15, no. 2 (2003): 221-231."
      ]
    },
    {
      "metadata": {
        "id": "OueLQM9xWZpF",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "\n",
        "\n",
        "# define function calculating derivatives\n",
        "def func_act(variables, t, S=lambda time:50):\n",
        "    [R, E] = variables\n",
        "    k_R = [0.22, 0.001]\n",
        "    k_E = [0.8, 0.01]\n",
        "    Km = [0.01, 0.01]\n",
        "    E_total = 0.5\n",
        "    dR = k_R[0] * (E_total - E + S(t)) - k_R[1] * R\n",
        "    dE = k_E[0] * (E_total - E) / (E_total - E + Km[0])\\\n",
        "       - k_E[1] * R * E / (E + Km[1])\n",
        "    return [dR, dE]\n",
        "\n",
        "\n",
        "# define gaussian function\n",
        "def gaussian(x, a, c, s):\n",
        "    return a * np.exp(-(x - c)**2 / (2 * s**2) )\n",
        "\n",
        "\n",
        "S = lambda time:   gaussian(time, .1, 10000, 1000) \\\n",
        "                 + gaussian(time, .2, 30000, 1000) \\\n",
        "                 + gaussian(time, .3, 50000, 1000) \\\n",
        "                 + gaussian(time, .4, 70000, 1000)\n",
        "initial_position = [0, 0]\n",
        "t = np.arange(0.0, 100000.0, 1)\n",
        "\n",
        "diff_eqs = lambda vars, time: func_act(vars, time, S)\n",
        "solution = odeint(diff_eqs, initial_position, t)\n",
        "result = np.transpose(solution)\n",
        "signal = [100 * S(_t) for _t in t]\n",
        "\n",
        "plots = [plt.plot(t, result[0])[0],plt.plot(t, signal)[0]]\n",
        "legends = [r'$[R]$', r'$100 \\times [S]$']\n",
        "plt.legend(plots, legends, loc=2)\n",
        "plt.xlabel(r'time, $t$, sec', fontsize=16)\n",
        "plt.ylabel(r'concentration, [X]',\\\n",
        "           fontsize=16)\n",
        "# save results\n",
        "plt.draw()\n",
        "plt.savefig('day_10_1.eps')\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "J7o5vc5bWixo",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **9.** Bio-switch with co-activation, bifurcation diagram\n",
        "\n",
        "For more details please see article:\n",
        "\n",
        "Tyson, John J., Katherine C. Chen, and Bela Novak. \"Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell.\" *Current opinion in cell biology* 15, no. 2 (2003): 221-231."
      ]
    },
    {
      "metadata": {
        "id": "KWHQ4wzCWqW0",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "\n",
        "# define function calculating derivatives\n",
        "def func_act(variables, t, S=0):\n",
        "    [R, E] = variables\n",
        "    k_R = [0.22, 0.001]\n",
        "    k_E = [0.8, 0.01]\n",
        "    Km = [0.01, 0.01]\n",
        "    E_total = 0.5\n",
        "    dR = k_R[0] * (E_total - E + S) - k_R[1] * R\n",
        "    dE = k_E[0] * (E_total - E) / (E_total - E + Km[0])\\\n",
        "       - k_E[1] * R * E / (E + Km[1])\n",
        "    return [dR, dE]\n",
        "\n",
        "initial_low = [0, 0]\n",
        "initial_high = [200, 0]\n",
        "t = np.arange(0.0, 100000.0, 1)\n",
        "S = np.arange(0, 0.5, 0.01)\n",
        "\n",
        "low_R = []\n",
        "high_R = []\n",
        "for _S in S:\n",
        "    diff_eqs = lambda vars, time: func_act(vars, time, _S)\n",
        "    solution = odeint(diff_eqs, initial_low, t)\n",
        "    low_R.append(solution[-1,0])\n",
        "    solution = odeint(diff_eqs, initial_high, t)\n",
        "    high_R.append(solution[-1,0])\n",
        "    \n",
        "Sthr = S[-1]\n",
        "for [_low_R,_high_R,_S] in np.transpose([low_R,high_R,S]):\n",
        "    if _high_R - _low_R < 0.001:\n",
        "        Sthr = _S\n",
        "        break    \n",
        "\n",
        "plots = [plt.plot(S, high_R)[0], plt.plot(S, low_R)[0]]\n",
        "legends = [r'$[R]_{min}^{stable}$',r'$[R]_{max}^{stable}$']\n",
        "plt.legend(plots, legends, loc=2)\n",
        "plt.xlabel(r'signal strength, [S]', fontsize=16)\n",
        "plt.ylabel(r'response, [R]', fontsize=16)\n",
        "plt.axvline(x=Sthr, color='k', ls='dashed')\n",
        "plt.text(Sthr, high_R[0], r'bifurcation', fontsize=20,\\\n",
        "         rotation=90)\n",
        "# save results\n",
        "plt.draw()\n",
        "plt.savefig('day_11_1.eps')\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "KIOTbG6vW6lu",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **10.** Bio-switch with co-repression\n",
        "\n",
        "For more details please see article:\n",
        "\n",
        "Tyson, John J., Katherine C. Chen, and Bela Novak. \"Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell.\" *Current opinion in cell biology* 15, no. 2 (2003): 221-231."
      ]
    },
    {
      "metadata": {
        "id": "tUyLK2y1XAGK",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "\n",
        "\n",
        "# define function calculating derivatives\n",
        "def func_act(variables, t, S=lambda time:50):\n",
        "    [R, E] = variables\n",
        "    k_R = [1, 1, 0.01, 0.01]\n",
        "    k_E = [0.8, 0.01]\n",
        "    Km = [0.01, 0.01]\n",
        "    E_total = 0.5\n",
        "    dR = k_R[0] + k_R[1]*S(t) - (k_R[2] + k_R[3]*E)*R\n",
        "    dE = k_E[0] * (E_total - E) / (E_total - E + Km[0])\\\n",
        "       - k_E[1] * R * E / (E + Km[1])\n",
        "    return [dR, dE]\n",
        "\n",
        "\n",
        "# define gaussian function\n",
        "def gaussian(x, a, c, s):\n",
        "    return a * np.exp(-(x - c)**2 / (2 * s**2) )\n",
        "\n",
        "\n",
        "S = lambda time:  gaussian(time, .01, 10000, 2000) \\\n",
        "                + gaussian(time, .02, 30000, 2000) \\\n",
        "                + gaussian(time, .03, 50000, 2000) \\\n",
        "                + gaussian(time, .05, 70000, 2000)\n",
        "initial_position = [0, 0.1]\n",
        "t = np.arange(0.0, 100000.0, 1)\n",
        "\n",
        "diff_eqs = lambda vars, time: func_act(vars, time, S)\n",
        "solution = odeint(diff_eqs, initial_position, t)\n",
        "result = np.transpose(solution)\n",
        "signal = [1000*S(_t) for _t in t]\n",
        "\n",
        "plots = [plt.plot(t, result[0])[0],plt.plot(t, signal)[0]]\n",
        "legends = [r'$[R]$', r'$1000 \\times [S]$']\n",
        "plt.legend(plots, legends, loc=2)\n",
        "plt.xlabel(r'time, $t$, sec', fontsize=16)\n",
        "plt.ylabel(r'concentration, [X]',\\\n",
        "           fontsize=16)\n",
        "# save results\n",
        "plt.draw();\n",
        "plt.savefig('day_12_1.eps');\n"
      ],
      "execution_count": 0,
      "outputs": []
    },
    {
      "metadata": {
        "id": "aPOMgTj3XApO",
        "colab_type": "text"
      },
      "cell_type": "markdown",
      "source": [
        "# **11.** Bio-switch with co-repression, bifurcation diagram\n",
        "\n",
        "For more details please see article:\n",
        "\n",
        "Tyson, John J., Katherine C. Chen, and Bela Novak. \"Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell.\" *Current opinion in cell biology* 15, no. 2 (2003): 221-231."
      ]
    },
    {
      "metadata": {
        "id": "mrItw5ufXK0G",
        "colab_type": "code",
        "colab": {}
      },
      "cell_type": "code",
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from scipy.integrate import odeint\n",
        "\n",
        "# define function calculating derivatives\n",
        "def func_act(variables, t, S=50):\n",
        "    [R, E] = variables\n",
        "    k_R = [1, 1, 0.01, 0.01]\n",
        "    k_E = [0.8, 0.01]\n",
        "    Km = [0.01, 0.01]\n",
        "    E_total = 0.5\n",
        "    dR = k_R[0] + k_R[1]*S - (k_R[2] + k_R[3]*E)*R\n",
        "    dE = k_E[0] * (E_total - E) / (E_total - E + Km[0])\\\n",
        "       - k_E[1] * R * E / (E + Km[1])\n",
        "    return [dR, dE]\n",
        "\n",
        "initial_low = [0, 0]\n",
        "initial_high = [200, 0]\n",
        "t = np.arange(0.0, 100000.0, 1)\n",
        "S = np.arange(0, 0.1, 0.01)\n",
        "\n",
        "low_R = []\n",
        "high_R = []\n",
        "for _S in S:\n",
        "    diff_eqs = lambda vars, time: func_act(vars, time, _S)\n",
        "    solution = odeint(diff_eqs, initial_low, t)\n",
        "    low_R.append(solution[-1,0])\n",
        "    solution = odeint(diff_eqs, initial_high, t)\n",
        "    high_R.append(solution[-1,0])\n",
        "\n",
        "Sthr = S[-1]\n",
        "for [_low_R,_high_R,_S] in np.transpose([low_R,high_R,S]):\n",
        "    if _high_R - _low_R < 0.001:\n",
        "        Sthr = _S\n",
        "        break\n",
        "\n",
        "plots = [plt.plot(S, high_R)[0], plt.plot(S, low_R)[0]]\n",
        "legends = [r'$[R]_{min}^{stable}$',r'$[R]_{max}^{stable}$']\n",
        "plt.legend(plots, legends, loc=2)\n",
        "plt.xlabel(r'signal strength, [S]', fontsize=16)\n",
        "plt.ylabel(r'response, [R]', fontsize=16)\n",
        "plt.axvline(x=Sthr, color='k', ls='dashed')\n",
        "plt.text(Sthr, high_R[0], r'bifurcation', fontsize=20,\\\n",
        "         rotation=90)\n",
        "# save results\n",
        "plt.draw();\n",
        "plt.savefig('day_13_1.eps');"
      ],
      "execution_count": 0,
      "outputs": []
    }
  ]
}