2020_Kinetic_modeling_in_Python_student.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "anaconda-cloud": {},
    "colab": {
      "name": "2020_Kinetic_modeling_in_Python_student.ipynb",
      "provenance": [],
      "collapsed_sections": [
        "W9cNb5nJch13"
      ]
    },
    "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.3"
    }
  },
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "IwQ2BMkechyR"
      },
      "source": [
        "<h1 align=\"center\"><font color='darkgreen'> Practical: kinetic modeling in Python </font></h1>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "_tVNTY8ZchyS"
      },
      "source": [
        "Oct 2020 \n",
        "    \n",
        "Teaching Assistants: <b>Douwe Grundel, Aleksandr Pogodaev</b>\n",
        "\n",
        "(Molecular Systems Biology group, RUG)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "jBEiHNXKchyT"
      },
      "source": [
        "**This practical has three goals:**\n",
        "\n",
        "* familiarizing with Python that offers a plethora of tools to process, analyse and plot scientific data;\n",
        "\n",
        "* familiarizing with the basics of kinetic modeling and numerical integration of systems of ordinary differential equations. \n",
        "\n",
        "* A practical exercise to get a feeling of the added value that kinetic models can provide in life science research\n",
        "\n",
        "\n",
        "First, we will recapitulate basic chemical kinetics theory (reactions of first and second order); we will analyze the equations by hand and solve them numerically with one of the Python modules.  \n",
        "\n",
        "\n",
        "Second, we will work on a small model of cell cycle control and a model of a flux-dependent oscillatory system. Lastly, we will analyse models of a bistable system and the blood coagulation pathway. These examples should provide you with an idea how deeper understanding about a biological system can be attained by dynamic modeling in comparison to the static depiction of biological networks classically presented in textbooks of cell biology, biochemistry and the like.\n",
        "\n",
        "\n",
        "**Grade:**\n",
        "For each exercise, it is indicated how many points can be earned for the correct answer. The total of **111** points can be obtained, so your grade for this practical will be: $\\frac{scored \\; points}{111} \\;* 10$.\n",
        "Please make all exercises **in this file!** At the end of the practical, send the complete file to d.a.j.grundel@rug.nl"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "JZJ4SIn2chyU"
      },
      "source": [
        "***Double-click*** here to enter the edit mode.\n",
        "\n",
        "*Name:*\n",
        "\n",
        "*Student number:*\n",
        "\n",
        "Use ***Shift+Enter*** to leave the edit mode."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "NhaTHC3NchyV"
      },
      "source": [
        "<h2 align=\"left\"><font color='navy'> A. Mass action kinetics [24 points]</font></h2>\n",
        "<h3 align=\"left\"><font color='navy'><i>Model 1. Irreversible first order reactions [8 points]</i></font></h3>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "bDKxFFUzchyW"
      },
      "source": [
        "First-order reactions are reactions in which one substrate $S$ is irreversibly transformed into product $P$. In generalized schematics, the reaction looks the following:\n",
        "\n",
        "<center>$S \\xrightarrow[]{r_1} P$</center>\n",
        " \n",
        "Prominent examples are decay reactions like the thermic decay of $N_2O_5$: \n",
        "\n",
        "<center>$ N_2O_5 \\xrightarrow[]{r_1} N_2O_4+\\frac{1}{2}O_2$</center>\n",
        " \n",
        "or radioactive decay (depicted here is the alpha-decay of Radium to Radon): \n",
        "\n",
        "<center>$ Ra \\xrightarrow[]{r_1} Rn+ 𝛼$</center>\n",
        " \n",
        "The rate of the reaction is the temporal change (*i.e.* the derivative of the concentration with respect to time) in substrate or product:\n",
        "\n",
        "<center>$\\frac{dS}{dt} = −r_1$</center>\n",
        " \n",
        "<center>$\\frac{dP}{dt} = r_1$</center>\n",
        " \n",
        "**Mass action** kinetics assumes that this rate is merely dependent on the concentration of substrate present: \n",
        "\n",
        "<center>$r_1=k_1\\cdot S$,</center>\n",
        " \n",
        "with $S$ representing the concentration of substrate and $k_1$ – the constant of proportionality which is usually called the rate constant.\n",
        "<br>\n",
        "<br>\n",
        "<br>\n",
        "<br>\n",
        "\n",
        "\n",
        "<font color='firebrick'><b>Task A1: Find the analytic solution for the differential equation that describes the change in $S$ [1 points].</b></font>\n",
        "\n",
        "*The goal is to derive the expression $S(t) = f(t)$, where $f$ is a function of time. Solve by separation of variables!*\n",
        "\n",
        "*You need the following integration rule:* "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "fLINykHDchyX"
      },
      "source": [
        "<center>$\\int \\frac{1}{x}dx = ln|x|+C$</center>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "E6DneTA2chyY"
      },
      "source": [
        "<font color='firebrick'><i>*Write your analytic derivation here using $\\LaTeX$ code.*</i></font> *If you double click on the integration formula given above you will see the $\\LaTeX$ code used to write it. From it, get the idea how to write your derivation. The powerful google may help you as well. You can return from the coding mode to the output by cellecting the cell of iterest and hitting Shift+Enter.*\n",
        "\n",
        "\n",
        "$...$\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "R7LtJgmRchyZ"
      },
      "source": [
        "The analytical solution gives the concentrations of $S$ and $P$ as a function of time. For a small system, the differential equations are relatively easy to solve. However, when a system involves more reactions it quickly gets very complicated. Therefore, it is often more practical to solve the system numerically. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "i9LsOwwBchya"
      },
      "source": [
        "<font color='firebrick'><b>Task A2: Integrate the differential equation numerically [0 points].</b></font>\n",
        "\n",
        "*Use the following values for the initial concentrations and the parameter: \n",
        "$S = 1\\ mM$, $P = 0\\ mM$, $k_1 = 0.5\\ s^{-1}$*\n",
        "\n",
        "\n",
        "**Solving a system numerically** works as follows: After you write the code and hit run, a software fills in the values for the parameters and the initial concentrations in the rate equations, and calculates the rates. Then it takes a small step in time ($\\delta t$), while assuming that the reaction rates stay constant during this time step. The concentrations at time $t + \\delta t$ are determined using these reaction rates. Subsequently, the concentrations at $t + \\delta t$ are used to calculate the reaction rates again. Then another time step is taken and the concentrations at time $t + 2\\delta t$ are determined. This process continues until the predetermined final time point is reached. \n",
        "\n",
        "\n",
        "The assumption that the reaction rates are constant during the interval $\\delta t$ is, of course, not true. The rate equations are functions of the concentrations of the reactants, and these concentrations change during the time step. Nonetheless, by taking a very small value for $\\delta t$, it is possible to get an accurate result. The Python module that we will use contains a number of algorithms (which we will not discuss here) that take into account that the reaction rates are not constant during a time step, thereby decreasing (but not completely eliminating) the error. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "2uAdXUBqchyb"
      },
      "source": [
        "In the following code it will be shown how to solve differential equations numerically using Python's module Tellurium.\n",
        "\n",
        "#### 0. Importing required modules ####"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "A0hADyAachyc"
      },
      "source": [
        "We start a script with importing the mathematical functions and classes that we will use. Python itself is just a regular programming language, and doesn't know about the Tellurium library until you ask to import it. This first statement is an easy way to get everything imported at once (to run a code in the cell of interest select it and hit Shift+Enter or press the Play button (triangle) on the left side of the cell), wait a bit until the cell is executed:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "cfFVT1KHchyd",
        "colab": {}
      },
      "source": [
        "!pip install tellurium &> /dev/null #installing the module into your local Colab environment\n",
        "import tellurium as te              #importing the module for this particular notebook"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "-T5tmGWrchyj"
      },
      "source": [
        "Now let's import some modules required for data processing and plotting:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "CiDNkYkbchyj",
        "colab": {}
      },
      "source": [
        "#Importing modules for imaging and data processing\n",
        "\n",
        "%matplotlib inline\n",
        "\n",
        "import matplotlib\n",
        "import numpy as np\n",
        "import matplotlib.pyplot as plt"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "v8Xcd7a3chym"
      },
      "source": [
        "*Python note: any line beginning with one or more **#** characters is a comment in Python. A comment may also be added to the end of a line, after which no more Python code can follow until the next line.*"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "JfxP4jng0I1-"
      },
      "source": [
        "*If Python returns you a large message with an error, first rule: **don't panic**. Mostly probably there is a typo in your code. Read the error message, it will point you to the place where the error appears. Then carefully check your code. If you are struggling for some time, then ask your fellow students or tutorial supervisors to help. Use of a comma instead of a dot or vice versa could often be the source of the problem.*"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Y0DbJv00chyn"
      },
      "source": [
        "**Generally, in order to perform a dynamic simulation one has to accomplish the following actions:**\n",
        "1. identify the variables and give their initial values (in context of biology, concentrations of compound are variables),\n",
        "2. define reactions,\n",
        "3. define their rate equations,\n",
        "4. set values for the kinetic parameters.\n",
        "\n",
        "In the Tellurium package a model is written in a text string, that contains the information mentioned in #1-4.\n",
        "Let's take a look at an example of a model string below."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "-whX9DihR5vY",
        "colab": {}
      },
      "source": [
        "  model_stringA2='''\n",
        "  # Simple model of a first order reaction\n",
        "  # Written by RUG Student, 2020\n",
        "  \n",
        "  # Identify variables and specify initial values:\n",
        "  S = 1;  P = 0;\n",
        "\n",
        "  # Reactions with rate laws:\n",
        "  R1: S -> P; k1*S;\n",
        " \n",
        "  # Kinetic parameters:\n",
        "  k1 = 1.0;\n",
        "'''"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "X3cGmkkfEv9w"
      },
      "source": [
        "From this string you can see, that \n",
        "- a model string in Tellurium should be specified within three quotations, i.e. **model string = ''' ... '''**\n",
        "- comments, that help to structure the string are marked with **#**. You may add text comments to the model to make it clearer, but also additinally to describe what you model, who are you, when you wrote the model, etc.\n",
        "- variables are specified together with initial conditions, i.e. **S=1**. Here values are separated by a semicolon **;**\n",
        "- reactions should follow this pattern **reaction name: reactants -> products; rate law;**\n",
        "- finally kinetic parameters (rate constants) are defined, i.e. **k1=1.0**, separated by a semicolon **;**\n",
        "\n",
        "**NB:** It's important to keep the units in the model consistent, i.e. don't use both µM and mM simultaneosly, only one type of units"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "6YtR_pv0chy8"
      },
      "source": [
        "#### 5. Passing the model_string to Tellurium ####"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "77EfKQJ60I2I"
      },
      "source": [
        "The next step is to pass the created string to Tellurium. Tellurium in turn will create a set of Ordinary Differential Equations (ODEs) automatically."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "pSu0BXvWchzG",
        "colab": {}
      },
      "source": [
        "modelA2 = te.loada(model_stringA2)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "HS8QUKxrA-zP"
      },
      "source": [
        "#### 6. Performing the simulation ####"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "MeaZ4lSlchzJ"
      },
      "source": [
        "We are now equipped to solve the ODE to obtain a trajectory from initial conditions. \n",
        "\n",
        "A good practice when using Tellurium is to reset the model before performing the simulation, using ```modelA2.resetAll()```. It ensures that the model is in the initial state.\n",
        "We simulate the model by calling ```modelA2.simulate(start=0, end=10, points=100) ```. Where ```start``` and ```end``` specify start end time of simulation respectively, ```points``` specify how many integration points will be between the start and end time. We store the result of simulation in the variable ```resultA2```.\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "c3F4AMHlchzK",
        "colab": {}
      },
      "source": [
        "modelA2.resetAll()\n",
        "resultA2=modelA2.simulate(start=0, end=10, points=100)    "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "ZmH0AZUOchzO"
      },
      "source": [
        "#### 7. Plotting the trajectories ####"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "K24mxtIO0I2Y"
      },
      "source": [
        "There are two ways to plot the data from the simulation. First is to use in-built feature of Tellurium like that:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "VeqqsCjrchzR",
        "colab": {}
      },
      "source": [
        "modelA2.plot(xtitle='Time, min',ytitle='Concentrations, mM')      #Plot the results"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "I0JsZ1G30I2d"
      },
      "source": [
        "The second option is to use a popular Python module called *matplotlib*. Here we take the simulated concentrations and time from the result obtained above. By calling ```resultA2['time']``` we get time points and by calling ```resultA2['[S]']``` we get calculated concentrations of $S$, the same principle works to get concentrations of $P$."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "rEdqkTvn0I2e",
        "colab": {}
      },
      "source": [
        "plt.plot(resultA2['time'],resultA2['[S]'],label='[S]') #Plot the line for S\n",
        "plt.plot(resultA2['time'],resultA2['[P]'],label='[P]') #Plot the line for P\n",
        "\n",
        "plt.legend()                                           #Call the legend\n",
        "plt.xlabel('Time, min')                                #Specify label for the x axis\n",
        "plt.ylabel('Concentration, mM')                        #Specify label for the y axis\n",
        "plt.show()                                             #Show the plot"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "qPCJvLJA0I2j"
      },
      "source": [
        "While you need to write more lines of code to plot with matplotlib, it offers more tunability than inbuilt plot() feature of Tellurium. Later in this tutorial you will see that we use matplotlib to plot additional more complicated graphs.\n",
        "\n",
        "**Important:** Always label axes on your plots and add legend! We will subtract points for not annotated plots."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "nzYT7UxXchzX"
      },
      "source": [
        "<font color='firebrick'><i>Question A2.1:Why does the rate at which $P$ is produced decrease over time? [1 point]</i></font> \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "hOWrf6suBwbN"
      },
      "source": [
        "Answer A2.1: "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "77AssvVh7oQI"
      },
      "source": [
        "### Changing parameters and initial conditions. ###"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "H62Ez5Fcc4QG"
      },
      "source": [
        "We have written ```modelA2``` and specified default initial conditions and parameters in the strings. To see the effect of different initial conditions and parameters you don't need to rewrite the whole model. Let's explore how we can do it.\n",
        "\n",
        "After each simulation, Tellurium stores the final concentration as initial conditions for the next simulation. If you want to perform a new simulation with your default initial conditons , you need to reset the model by calling:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "6hcbMKrYdZyk",
        "colab": {}
      },
      "source": [
        "modelA2.resetAll()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "ucW8Gd2icVvl"
      },
      "source": [
        "To change reaction rates or starting concentrations we don't need to rewrite the whole model. Instead we can change them calling them from the model:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "-4I7hBd2crMx",
        "colab": {}
      },
      "source": [
        "modelA2.k1=0.5  #Change kinetic parameter k1 from k1=1.0 to k1=0.5\n",
        "modelA2.S=2     #Change the initial concentration of S from S=1 to S=2"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Y2U0a32WeOVn"
      },
      "source": [
        "And then we simulate the model with updated reaction rate k1 and initial concentration S1. What has changed from the previous simulation on the graph?"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "FN-FxvK5eEyC",
        "colab": {}
      },
      "source": [
        "modelA2.simulate(start=0, end=10, points=100)                     #Simulate the model\n",
        "modelA2.plot(xtitle='Time, min',ytitle='Concentrations, mM')      #Plot the results"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "k0g6kgvm8JCn"
      },
      "source": [
        "<font color='firebrick'><i>Question A2.2: In the previous graph  around $1.38$ minutes the concentration curve of $S$ intersects that of $P$. If the value of $k_1$ is decreased to $0.05$, what would happen to the time of intersection? Test your hypothesis with the model!  Tip: Increase end time of simulation if needed.\n",
        "[2 points]</i></font>  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "dN5S7SuN-66X"
      },
      "source": [
        "Answer A2.2: \n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "yu78W4spchzY",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE MODEL IN THE QUESTION A2.2 HERE.\n",
        "#FOLLOW THE SEQUENCE RESET - NEW PARAMETERS - SIMULATE - PLOT\n",
        "\n",
        "  "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "0jAmUHzQchzi"
      },
      "source": [
        "<font color='firebrick'><i>Question A2.3:The final concentration of $P$ is $1$. What would the final concentration of $P$ be if the initial concentrations of $S$ and $P$ were $2$ and $3$, respectively? [2 points]</i></font>  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "LW9HwDRAGS4M"
      },
      "source": [
        "Answer 2.3: "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "zmrNrmUrchzj",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE MODEL IN THE QUESTION A2.3 HERE:\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "CE0h1qkhchzq"
      },
      "source": [
        "<h3 align=\"left\"><font color='navy'><i>Model 2. Reversible first order reactions [9 points]</i></font></h3>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "VSqNQ-Eqchzr"
      },
      "source": [
        "Most chemical (and biochemical) reactions are not irreversible. We will therefore expand the model with a second reaction, the reverse of $r_{r}$, where product $P$ turns back into substrate $S$. The reverse reaction follows mass action kinetics too. In schematics, this looks as follows:\n",
        "\n",
        "<center>$ S \\underset{k_{r}}{\\stackrel{k_f}{\\rightleftharpoons}} P $</center>\n",
        " \n",
        "<font color='firebrick'><b>Task A3: Write down the set of differential equations describing this reaction pair. Solve it numerically [3 points].</b></font>\n",
        "\n",
        "\n",
        "*Follow the same steps as you did for **Model 1** to incorporate the second reaction into your model. Take a value of $0.5\\ s^{-1}$ for $k_f$, of $0.2\\ s^{-1}$ for $k_{r}$, and of $1\\ mM$ for $S$ at the start of the reaction.*\n",
        "\n",
        "There are two ways to write a reversible reaction in Tellurium. You can use either of them.\n",
        "```\n",
        "First way                     Second way:\n",
        "R1: S -> P; k_f*S             R1: S -> P; k_f*S - k_r*P  #Note that -> arrow remains the same,\n",
        "R2: P -> S; k_r*P                                        #even though we describe a reverisble reaction\n",
        "``` "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "VE55nWWvchzs",
        "colab": {}
      },
      "source": [
        "#WRITE THE MODEL STRING HERE AND PASS IT TO TELLURIUM\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "reWpgxHBchzw",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE RESULTS HERE\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "3TPuJVdBchz4"
      },
      "source": [
        "<font color='firebrick'><i>Question A3.1: In your own words, explain why the concentration of $S$ does not approach $0$ [3 points].</i></font>\n",
        "<br>\n",
        "<br>\n",
        "<font color='firebrick'><i>Question A3.2: Report the steady state concentration of $S$ and $P$ in the simulation results. How do you know the system reached steady state? [3 points].</i></font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "hPgxQm-tchz5"
      },
      "source": [
        "<h3 align=\"left\"><font color='navy'><i>Model 3. Reversible second order reactions [9 points].</i></font></h3>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "MiU8rM_Ochz6"
      },
      "source": [
        "In second order reactions two substrates react to one product. The model can be schematically represented as follows: \n",
        "\n",
        "<center>$ S_1 + S_2 \\underset{k_{r}}{\\stackrel{k_f}{\\rightleftharpoons}} P $</center>\n",
        " \n",
        "<font color='firebrick'><b>Task A4: Write down the set of differential equations describing this reaction pair. Solve it numerically [3 points].</b></font>\n",
        "\n",
        "\n",
        "*Follow the same steps as you did for **Model 1** and **2** to incorporate the second substrate into your model. Take a value of $0.5\\ mM^{-1} s^{-1}$ for $k_f$ and $0.2\\ s^{-1}$ for $k_{r}$. Set the initial amount of $S_1$ to be $1\\ mM$, of  $S_2$ to be $0.8\\ mM$ and of $P$ to be $0\\ mM$.*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "toB7EGH5chz7",
        "colab": {}
      },
      "source": [
        "#WRITE THE MODEL STRING HERE AND PASS IT TO TELLURIUM\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "OSclx9pHchz_",
        "colab": {}
      },
      "source": [
        "#PLOT THE RESULTS HERE\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "_qWgiXyZch0G"
      },
      "source": [
        "<font color='firebrick'><i>Question A4.1: Find the steady state concentrations of $S_1$, $S_2$ and $P$ in your simulation [3 points].</i></font>\n",
        "<br>\n",
        "<br>\n",
        "<font color='firebrick'><i>Question A4.2: Change the initial concentration of $S_2$ to $0.1$. Rerun the simulation and find the steady state concentrations of $S_1$, $S_2$ and $P$. Also try an initial $S_2$ concentration of $5$. In your own words, describe the differences and explain the cause [3 points].</i></font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "152izmZ9K87h"
      },
      "source": [
        "Answer 4.1:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "r--hYePxLPu5"
      },
      "source": [
        "Answer 4.2:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "4uSJYjqach0J",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE MODEL IN THE QUESTION A4.2 HERE:\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "vq1j5IFPch0S",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE MODEL IN THE QUESTION A4.2 HERE:\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "K8pJXgoich0Y"
      },
      "source": [
        "<h2 align=\"left\"><font color='navy'> B. Cell cycle control [27 points]</font></h2>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "2kRY2g4cch0a"
      },
      "source": [
        "The cell cycle is the series of events that lead to division and duplication of a cell. The molecular events that control the cell cycle are ordered and directional; that is, each process occurs in a sequential fashion and it is impossible to reverse the cycle. Therefore, the regulatory components need to oscillate in time. \n",
        "\n",
        "\n",
        "We assume the cell cycle to be driven by a protein circuit centered on the cyclin-dependent protein kinase (CDK1) and the anaphase-promoting complex (APC). CDK1 is activated by the rapid, high-affinity binding of cyclin, which itself is being synthesized at a constant rate. Activated CDK1 phosphorylates a variety of target substrates driving the cell into mitosis. It also activates the APC which drives the cell again out of mitosis and therefore needs to lag behind CDK1. APC will inactivate CDK1. The cycling appearance of CDK1 and APC is depicted in **Figure 1**.  "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "lXlXlM6dch0b"
      },
      "source": [
        "<img src=\"https://github.com/vakiltakhaveev/misc/blob/master/cell_cycle_osc.png?raw=true\" align=\"center\" width=\"300\">\n",
        "\n",
        "**Figure 1**: CDK 1 drives the cell into mitosis, APC lagging behind drives it out again.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "NBJzpH15ch0c"
      },
      "source": [
        "Here, we try to understand what is necessary for this cycling, oscillatory appearance of the different components. Therefore, we are going to analyze three possible models of rising complexity of cell cycle control. They are depicted in **Figure 2**."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "E6EbgJllch0d"
      },
      "source": [
        "<img src=\"https://github.com/vakiltakhaveev/misc/blob/master/cell_cycle_models.png?raw=true\" align=\"center\" width=\"600\">\n",
        "\n",
        "**Figure 2**: Different negative feedback loops of rising complexity. **A.** CDK1 exerts negative feedback on itself. **B.** CDK1 activates APC which leads to its own inactivation. **C.** The polo-like kinase 1 has been added to the regulatory network. It gets activated by CDK1 and then activates APC. In turn, APC inactivates CDK1. For all three it is assumed that CDK1 itself is activated by the rapid, high-affinity binding of cyclin, which is being synthesized at a constant rate. "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "frRex4Vech0e"
      },
      "source": [
        "<h3 align=\"left\"><font color='navy'><i>Model 1. Cell cycle control mediated only by CDK1 [4 points]</i></font></h3>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "4aFWC6pWch0f"
      },
      "source": [
        "We first analyze the simple one-component model depicted in **Figure 2A**. It can be described the following: CDK1 is activated by the binding of cyclin, which is present at a constant concentration. Because the cyclin concentration is constant, we assume that the activation rate of CDK1 is also constant:\n",
        "\n",
        "<center>$ r_{activation} = k_1$</center>\n",
        " \n",
        "The inactivation of CDK1 could be described as a mass action kinetic of second order depending on CDK1 and APC: \n",
        "\n",
        "<center>$ r_{inactivation} = k_2 \\cdot CDK1 \\cdot APC$</center>\n",
        " \n",
        "However, here we want to construct a very simple model based only on CDK1 and, hence, we want to remove APC from the rate equation. Therefore, we assume that APC is regulated so fast, that we can view it as a direct function of CDK1. We will assume that this regulation is ultrasensitive and describe it by a Hill function.\n",
        "\n",
        "<center>$ r_{inactivation} = k_2 \\cdot CDK1 \\cdot \\frac{CDK1^n}{K^n+CDK1^n}$</center>\n",
        " "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "VYuHJwLLch0h"
      },
      "source": [
        "<font color='firebrick'><b>Task B1: Numerically integrate the equation over time [4 points].</b></font>\n",
        "\n",
        "\n",
        "*Use the following initial concentration and parameters:*\n",
        "\n",
        "\n",
        "*$CDK1 = 0$, $k_1 =  0.1$, $k_2 = 1$, $K = 0.5$, $n =8$*"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "b4fYWCAlch0i",
        "colab": {}
      },
      "source": [
        "#WRITE THE MODEL STRING HERE AND PASS IT TO TELLURIUM\n",
        "#USE ^ for the power, e.g.: CDK1^n\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "SjHVqB-gch0l",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE RESULTS HERE\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "5SVpJlM5ch0o"
      },
      "source": [
        "<h3 align=\"left\"><font color='navy'><i>Model 2. Cell cycle control mediated by CDK1 and APC [11 points]</i></font></h3>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Sgisimn_ch0o"
      },
      "source": [
        "In the next step, we add the second ODE, which describes that APC regulation is not instantaneous (**Figure 2B**).\n",
        "<br>\n",
        "<br>\n",
        "\n",
        "\n",
        "<font color='firebrick'><b>Task B2.1: Write a model consisting of the two components, CDK1 and APC, changing over time and simulate it [3 points]. </b></font>\n",
        "\n",
        "\n",
        "*Give the initial conditions of CDK1 and APC: both start with an active concentration of zero. Assume again that CDK1 is activated at a constant rate. For the inactivation of CDK1 assume that APC acts in an ultrasensitive process.\n",
        "APC gets activated by the ultrasensitive binding of CDK1. Assuming that the total concentration of APC is constant, we take the concentration of inactive APC to be $1-APC_{active}$. Put together, the activation for APC looks the following:*\n",
        "\n",
        "\n",
        "*<center>$ r_{activation\\_APC} = k \\cdot (1-APC) \\cdot \\frac{CDK1^n}{K^n+CDK1^n} $</center>*\n",
        "\n",
        "\n",
        "*Assume that APC gets inactivated by a process following mass action kinetics of 1st order.*\n",
        "\n",
        "\n",
        "*Use the following parameters for the simulation:*\n",
        "\n",
        "\n",
        "*$k_1 = 0.1$ (activation of CDK1),*\n",
        "\n",
        "*$k_2 = 3$ (inactivation of CDK1),*\n",
        "\n",
        "*$k_3 = 3$ (activation of APC),*\n",
        "\n",
        "*$k_4 = 1$ (inactivation of APC),*\n",
        "\n",
        "*$n_1 = 8$ (Hill coefficient for the action of APC onto CDK leading to CDK’s inactivation),*\n",
        "\n",
        "*$n_2 = 8$ (Hill coefficient for the action of CDK onto APC leading to APC’s activation),*\n",
        "\n",
        "*both $K = 0.5$.*\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "3F6J6ibEch0q",
        "colab": {}
      },
      "source": [
        "#WRITE THE MODEL STRING HERE AND PASS IT TO TELLURIUM\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "Zu8QQG1Sch0t",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE RESULTS HERE\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "ek-HZAGYch0v"
      },
      "source": [
        "<font color='firebrick'><b>Task B2.2: Use matplotlib (see above) to depict the system in a phase space. That is, plot CDK1 on the x-axis and APC on the y-axis. The phase space represents all possible states of a system. The system’s evolving state over time traces a path (a phase space trajectory).[2 points]</b></font>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "JQ-bnxMNch0y",
        "colab": {}
      },
      "source": [
        "#PLOT THE PHASE SPACE HERE\n",
        "                                           "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "H4UjVKpMXsw8"
      },
      "source": [
        " <font color='firebrick'><i> Question B2.1. What is the final state of the system in the phase space plot? Can you find the same in the time course plot? [2 points]</i></font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Buady2IwUhCn"
      },
      "source": [
        "Answer B2.1: \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "8aFc_8W9ch01"
      },
      "source": [
        "<font color='firebrick'><b>Task B2.3: Now let's check if differing starting conditions lead you to a different steady-state. Create a phase space figure (again with matplotlib) containing trajectories starting from several different inial conditions (see cell below)  .[2 points]</b></font>\n",
        "\n",
        "   "
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "AYjXQt5gch02",
        "colab": {}
      },
      "source": [
        "#let's make a list of pairs of initial conditions (CDK1, APC)\n",
        "pairs_of_initial_conditions = [(0, 0), (0.1, 0.35), (0.3, 0.45), (0.55, 0.05), (0.5, 0.4)]\n",
        "\n",
        "#no we will simulate the model several times; each time with a different pair of initial conditions\n",
        "#for that we will iterate over all pairs in the list of initial conditions that we created above\n",
        "\n",
        "for pair in pairs_of_initial_conditions:\n",
        "    modelB2.resetAll() # if you used a different name for your model, you will need to replace 'modelB2' with that name\n",
        "\n",
        "    modelB2.CDK1 = pair[0] # pair[0] is the first element of the pair\n",
        "    modelB2.APC = pair[1]  # pair[1] is the second element of the pair\n",
        "    \n",
        "    # FINISH THE CODE AND MAKE SURE EACH INITIAL CONDITION IN pairs_of_initial_conditions IS PLOTTED"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "5l6l1rftch06"
      },
      "source": [
        "<font color='firebrick'><i>Question B2.2: How could you interpret this figure? [2 points]</i></font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "aEvDc2v1Z_KK"
      },
      "source": [
        "Answer B2.2:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "f7Gj5Bqxch07"
      },
      "source": [
        "<h3 align=\"left\"><font color='navy'><i>Model 3. Cell cycle control model with three components [12 points]</i></font></h3>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "vB4MIx2Gch07"
      },
      "source": [
        "<font color='firebrick'><b>Task B3.1: Add a third reactant to the model (Plk1, look at **Figure 2C**). Suppose Plk1 to act as an intermediate between CPK1 and APC [3 points]. </b></font>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "8iCSmzDGch08"
      },
      "source": [
        "*Suppose that CDK1 gets activated and inactivated as before. Plk1 is activated by CDK1 (assume ultrasensitivity; concentration of inactive Plk1 is given by (1-Plk1); inactivation follows mass action kinetics). APC is activated by Plk1 (assume ultrasensitivity); the inactivation follows mass action kinetics. Increase integration time to 20.*\n",
        "\n",
        "\n",
        "\n",
        "*Simulate these ODEs. Use the following parameters: $k_1 = 0.1$ (activation of CDK1), $k_2 = 3$ (inactivation of CDK1), $k_3 = 3$ (activation of Plk1), $k_4 = 1$ (inactivation of Plk1), $k_5 = 3$ (activation of APC), $k_6 = 1$ (inactivation of APC), assume all Hill coefficients to be $8$ and all $K$ to be $0.5$.*\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "JiHQ_hjzch09",
        "colab": {}
      },
      "source": [
        "#WRITE THE MODEL STRING HERE AND PASS IT TO TELLURIUM\n",
        "\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "LnZlOCs9ch1A",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE RESULTS HERE. SIMULATE FROM t=0 to t=20\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "OW30yTSDch1C"
      },
      "source": [
        "<font color='firebrick'><b>Task B3.2: Plot the phase space of the model. The model contains 3 variables, therefore the phase space will be 3-dimensional. [1 points]. </b></font>\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "gymB1xxcch1E",
        "colab": {}
      },
      "source": [
        "#importing a module need for 3D plotting \n",
        "from mpl_toolkits.mplot3d import Axes3D\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "SG8wlSCxch1G",
        "colab": {}
      },
      "source": [
        "#plotting the phase space\n",
        "fig = plt.figure()\n",
        "ax = Axes3D(fig)\n",
        "\n",
        "#Guess how to plot three variables! (you need to add one line of code to plot the line)\n",
        "ax.set_xlabel('CDK')\n",
        "ax.set_ylabel('APC')\n",
        "ax.set_zlabel('Plk1')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "-d_9Koj-ch1I"
      },
      "source": [
        "<font color='firebrick'><b>Task B3.3: In the cells below you will need to modify your model, but keep the (plotted graphs of the) original, unmodified model above. \n",
        "\n",
        "Below we will no longer assume ultrasensitive responses, but Michaelis-Menten instead. Modify the model accordingly (you can do so by simply changing the Hill coefficients!). [3 points]</b></font>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "AHTXmm0lch1J",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE RESULTS HERE\n",
        "#SIMULATE FROM t=0 to t=20\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "N8LXylfgch1N"
      },
      "source": [
        "<font color='firebrick'><i>Question B3: Compare your results between the three different models created in section B. In which case(s) did the system oscillate? What was the effect of removing the ultrasensitivity in the modified version of the third model? Given this comparison, what are some characteristics are necessary for an oscillatory system [5 points]. </i></font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Pc026e2DdmcK"
      },
      "source": [
        "Answer B3:\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "PG8Aamzdch1O"
      },
      "source": [
        "<h2 align=\"left\"><font color='navy'> C. Flux-dependent oscillatory circuit [30 points]</font></h2>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "WxcOOvb3ch1P"
      },
      "source": [
        "The following assignment is based on a paper that was published several years ago. The authors described the construction of a circuit that is capable of oscillating, depending on the flux through glycolysis. The oscillatory circuit was constructed by placing two genes of *E. coli*’s acetate pathway under control of acetylphosphate, an intermediate in this pathway. In order to characterize this circuit, the authors developed a kinetic model, which we will reproduce here with some small alterations."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "vW4PJLeAch1Q"
      },
      "source": [
        "<font color='firebrick'><b>Task C1: Read the following text and write down the differential equations that describe the system.</b></font>\n",
        "\n",
        "The system (see **Figure 3**) is designed around acetate metabolism in *E. coli*. $AcetylCoA$ is produced in glycolysis with rate $r_{gly}$, which we consider to be constant, at a value of $0.05\\ mM\\ min^{-1}$. It is consumed in a number of reactions. For simplicity, all $AcetylCoA$ consuming reactions (except for one, see below) are modeled as one reaction, which we assume to follow mass action kinetics, with a rate constant $k$ of $10\\ min^{-1}$.\n",
        "\n",
        "\n",
        "In the remaining reaction, $AcetylCoA$ is also converted into $AcetylPhosphate$ by the enzyme phosphate acetyltransferase ($Pta$). This enzyme follows Michaelis-Menten kinetics (at an enzyme concentration of $1\\ \\mu M$: $k_M = 0.06\\ mM$, $v_{max} = 80\\ mM\\ min^{-1}$). In the next step in the pathway, $AcetylPhosphate$ is reversibly converted into $Acetate$. Both the forward and reverse reactions follow mass action kinetics ($k_{forward} = 1\\ min^{-1}$; $k_{reverse} = 1\\ min^{-1}$).\n",
        "\n",
        "\n",
        "$Acetate$ is the substrate of the enzyme acetylCoA synthetase ($Acs$), which exhibits Michaelis-Menten kinetics (at an enzyme concentration of $1\\ \\mu M$: $k_M = 0.1\\ mM$, $v_{max} = 0.8\\ mM\\ min^{-1}$). $Acs$ converts acetate into $AcetylCoA$. $Acetate$ is also secreted by the cell with mass action kinetics ($k = 0.005\\ min^{-1}$).\n",
        "\n",
        "\n",
        "To achieve oscillations in the system, the expression of the genes encoding $Pta$ and $Acs$ has been placed under control of $AcetylPhosphate$. The expression of $Acs$ is increased in the presence of $AcetylPhosphate$, which we will model with the Hill equation ($v_{max} = 2\\ \\mu M\\ min^{-1}$, $K_H = 10\\ mM$, $n = 2$). The expression of $Pta$ is inhibited indirectly: first, $AcetylPhosphate$ stimulates the expression of the Lac repressor $LacI$, which is also modeled with the Hill equation ($v_{max} = 0.1 \\mu M\\ min^{-1}$, $K_H = 10\\ mM$, $n = 2$). In turn, $LacI$ represses the expression of $Pta$, because it was placed under control of a Lac operator. For this repression, use the following rate equation:\n",
        "\n",
        "\n",
        "<center>$r=v_{max} \\cdot(1- \\frac{LacI^n}{LacI^n + K_H^n})$,</center>\n",
        "\n",
        "where $r$ is the rate of $Pta$ synthesis, $v_{max} = 2 \\mu M\\ min^{-1}$, $K_H = 0.001\\ \\mu M$ and $n = 2$.\n",
        "The proteins ($Pta$, $Acs$ and $LacI$) are degraded according to mass action kinetics ($k = 0.06 min^{-1}$).\n",
        "\n",
        "\n",
        "\n",
        "<img src=\"https://github.com/vakiltakhaveev/misc/blob/master/flux_syst.png?raw=true\" align=\"center\" width=\"400\">\n",
        "<b>Figure 3:</b> A schematic representation of the system taken from the original publication. The production of acetyl-CoA from sugars and fatty acids will be modeled with a single rate only. The ATP and phosphate molecules are not included in the model. The reactions at the bottom of the figure will be simplified in our model, such that acetate ($OAc^-$) is directly exported, without being converted into acetic acid ($HOAc$) first.\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "rKQIlsQpch1R"
      },
      "source": [
        "<font color='firebrick'><b>Task C2: Simulate the model and plot the results [10 points].</b></font>\n",
        "\n",
        "\n",
        "*Use 0 as the initial concentrations. Integrate for at least 1000, because it will take the system some time to reach a steady-state.*\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "C2XCW3ftch1S",
        "colab": {}
      },
      "source": [
        "#WRITE THE MODEL HERE AND LOAD IT TO TELLURIUM\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "qUW5y6oBEwFH",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE MODEL HERE\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "D1rNDg0fch1U",
        "colab": {}
      },
      "source": [
        "#REPLACE NAMES OF THE VARIABLES USED HERE TO ONES YOU USED IN THE SIMULATION\n",
        "\n",
        "plt.figure(1, (8, 8))\n",
        "\n",
        "ax = plt.subplot(2, 2, 1)\n",
        "plt.plot(resultC2['time'], resultC2['[AcP]'],label='AcP') # Acetyl Phosphate\n",
        "plt.plot(resultC2['time'], resultC2['[HOAc]'], label='HOAc') # Acetate\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 2)\n",
        "plt.plot(resultC2['time'], resultC2['[AcCoA]'],label='AcCoA',color=\"red\") # acetyl CoA\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 3)\n",
        "plt.plot(resultC2['time'], resultC2['[Pta]'],label='Pta',color=\"purple\") # Pta\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 4)\n",
        "plt.plot(resultC2['time'], resultC2['[Acs]'],label='Acs',color=\"violet\") # Acs\n",
        "plt.plot(resultC2['time'], resultC2['[LacI]'], label='LacI',color=\"darkred\") # LacI\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "plt.tight_layout()\n",
        "plt.show()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "qJ1zX0u5ch1X"
      },
      "source": [
        "<font color='firebrick'><b> Task C3: Run two more simulations to test the effect of the glycolytic flux on the occurrence of oscillations in the circuit. Use $r_{gly} = 0.01$ and $0.5$ [4 points]. </b></font>\n",
        "<br>\n",
        "<br>\n",
        "\n",
        "<font color='firebrick'><i>Question C3.1: You should see completely different behaviors of the system. Depending on the glycolytic flux the system either reaches a steady-state, or it keeps oscillating. Can you explain why this is the case? [6 points]</i></font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "n8GEPf7Nilee"
      },
      "source": [
        "Answer C3.1:\n",
        "\n",
        "\n",
        "\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "I9lAyuplch1X",
        "colab": {}
      },
      "source": [
        "#SIMULATE THE MODEL HERE\n",
        "# case r_gly = 0.01\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "AQAjOCRMch1a",
        "colab": {}
      },
      "source": [
        "#REPLACE NAMES OF THE VARIABLES USED HERE TO ONES YOU USED IN THE SIMULATION\n",
        "\n",
        "plt.figure(1, (8, 8))\n",
        "\n",
        "ax = plt.subplot(2, 2, 1)\n",
        "plt.plot(resultC3['time'], resultC3['[AcP]'],label='AcP')\n",
        "plt.plot(resultC3['time'], resultC3['[HOAc]'], label='HOAc')\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 2)\n",
        "plt.plot(resultC3['time'], resultC3['[AcCoA]'],label='AcCoA',color=\"red\")\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 3)\n",
        "plt.plot(resultC3['time'], resultC3['[Pta]'],label='Pta',color=\"purple\")\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 4)\n",
        "plt.plot(resultC3['time'], resultC3['[Acs]'],label='Acs',color=\"violet\")\n",
        "plt.plot(resultC3['time'], resultC3['[LacI]'], label='LacI',color=\"darkred\")\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "plt.tight_layout()\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "Q8Wd7eKWch1d",
        "colab": {}
      },
      "source": [
        "#SIMULATE THE MODEL HERE\n",
        "### the case of r_gly = 0.5\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "L1oCEo9nch1e",
        "colab": {}
      },
      "source": [
        "#REPLACE NAMES OF THE VARIABLES USED HERE TO ONES YOU USED IN THE SIMULATION\n",
        "\n",
        "plt.figure(1, (8, 8))\n",
        "\n",
        "ax = plt.subplot(2, 2, 1)\n",
        "plt.plot(resultC3['time'], resultC3['[AcP]'],label='AcP')\n",
        "plt.plot(resultC3['time'], resultC3['[HOAc]'], label='HOAc')\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 2)\n",
        "plt.plot(resultC3['time'], resultC3['[AcCoA]'],label='AcCoA',color=\"red\")\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 3)\n",
        "plt.plot(resultC3['time'], resultC3['[Pta]'],label='Pta',color=\"purple\")\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 4)\n",
        "plt.plot(resultC3['time'], resultC3['[Acs]'],label='Acs',color=\"violet\")\n",
        "plt.plot(resultC3['time'], resultC3['[LacI]'], label='LacI',color=\"darkred\")\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "plt.tight_layout()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "OAafn575ch1j"
      },
      "source": [
        "<font color='firebrick'><b> Task C4: When the researchers who introduced this system in *E. coli* added acetate to the growth medium, they found that the oscillations disappeared (the cells carried a fluorescent reporter so the oscillations could be observed). In your model, include the uptake of extracellular acetate [6 points].</b></font>\n",
        "<br>\n",
        "<br>\n",
        "*For the uptake reaction assume mass action kinetics, with $k = 0.005\\ min^{-1}$. Consider the extracellular acetate concentration to be constant, at $30\\ mM$. Use $0.5\\ mM\\ min^{-1}$ for the glycolytic flux. Integrate for a longer time period to see the effect over a longer time period.*\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "RFgDXwBich1k",
        "colab": {}
      },
      "source": [
        "#WRITE THE MODEL STRING HERE AND LOAD IT TO TELLURIUM\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "6zbYJCl6P2Mf",
        "colab": {}
      },
      "source": [
        "#SIMULATE THE MODEL HERE\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "-Tu0xMxTch1m",
        "colab": {}
      },
      "source": [
        "#REPLACE NAMES OF THE VARIABLES USED HERE TO ONES YOU USED IN THE SIMULATION\n",
        "\n",
        "plt.figure(1, (8, 8))\n",
        "\n",
        "ax = plt.subplot(2, 2, 1)\n",
        "plt.plot(resultC4['time'], resultC4['[AcP]'],label='AcP')\n",
        "plt.plot(resultC4['time'], resultC4['[HOAc]'], label='HOAc')\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 2)\n",
        "plt.plot(resultC4['time'], resultC4['[AcCoA]'],label='AcCoA',color=\"red\")\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 3)\n",
        "plt.plot(resultC4['time'], resultC4['[Pta]'],label='Pta',color=\"purple\")\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 4)\n",
        "plt.plot(resultC4['time'], resultC4['[Acs]'],label='Acs',color=\"violet\")\n",
        "plt.plot(resultC4['time'], resultC4['[LacI]'], label='LacI',color=\"darkred\")\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "plt.tight_layout()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "FUDwzb_-ch1t"
      },
      "source": [
        "<font color='firebrick'><i> Question C5: \tExplain why uptake of acetate suppresses oscillations in the system. Does an extracellular concentration of 20 mM also result in a stable steady-state? (Why?) [4 points]</i></font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "cv02k5R_pN4K"
      },
      "source": [
        "Answer C5:\n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "RkzUFMdach1u",
        "colab": {}
      },
      "source": [
        "#SIMULATE THE MODEL HERE\n",
        "# 20 mM case\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "oxZApfkYch1y",
        "colab": {}
      },
      "source": [
        "#REPLACE NAMES OF THE VARIABLES USED HERE TO ONES YOU USED IN THE SIMULATION\n",
        "\n",
        "plt.figure(1, (8, 8))\n",
        "\n",
        "ax = plt.subplot(2, 2, 1)\n",
        "plt.plot(resultC5['time'], resultC5['[AcP]'],label='AcP')\n",
        "plt.plot(resultC5['time'], resultC5['[HOAc]'], label='HOAc')\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 2)\n",
        "plt.plot(resultC5['time'], resultC5['[AcCoA]'],label='AcCoA',color=\"red\")\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 3)\n",
        "plt.plot(resultC5['time'], resultC5['[Pta]'],label='Pta',color=\"purple\")\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "ax = plt.subplot(2, 2, 4)\n",
        "plt.plot(resultC5['time'], resultC5['[Acs]'],label='Acs',color=\"violet\")\n",
        "plt.plot(resultC5['time'], resultC5['[LacI]'], label='LacI',color=\"darkred\")\n",
        "\n",
        "plt.legend()\n",
        "plt.xlabel('Time, min')\n",
        "plt.ylabel('Concentration, mM')\n",
        "\n",
        "plt.tight_layout()"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "zcc7ujgbLsYy"
      },
      "source": [
        "<h2 align=\"left\"><font color='navy'> D. Analysing a bistable system [19 points]</font></h2> "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "S_GIG_XmLsYy"
      },
      "source": [
        "In this task we will consider a biochemical system where molecule P is produced and broken down. In order to keep things simple, we will not consider which molecules are required for the synthesis of P, nor will we consider the products produced when P is degraded. Instead, we will assume that P is generated (out of nothing) with rate $r_{form} = 3 + \\frac{20}{1+2^{-3(P-5)}}$. P is broken down with rate $r_{deg} = 2.5 P$. Both rates have the unit $mM \\;min^{-1}$."
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "pBJv6EJXLsYz"
      },
      "source": [
        "<font color='firebrick'><b>Task D1a: Simulate and plot the system. Also determine the steady state concentration of P.</b> [4 points] </font><br> Use the initial concentration of P = 3.0 mM."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "PHiLPhYeLsYz",
        "colab": {}
      },
      "source": [
        "# DEFINE MODEL AND SET INITIAL CONDITIONS\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "HI8mbmsFLsY1",
        "colab": {}
      },
      "source": [
        "# SIMULATE MODEL AND PLOT RESULS\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "_lvZjT3lLsZA"
      },
      "source": [
        "<font color='firebrick'>Please provide the steady state value of P:\n",
        "\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "zVCBOIISLsZA"
      },
      "source": [
        "<font color='firebrick'><b>Task D1b: Repeat task D1a with different initial condition [1 point].</b></font><br>\n",
        "Use the initial concentration of P = 7.0 mM"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "LDIgSmHjLsZB",
        "colab": {}
      },
      "source": [
        "# CHANGE INITIAL CONDITIONS\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "F8vWMTDTLsZC",
        "colab": {}
      },
      "source": [
        "# RUN SIMULATION AND PLOT RESULTS\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "T5l_F2UWLsZI"
      },
      "source": [
        "<font color='firebrick'>Please provide the steady state value of P:</b></font> "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "WhP52jl5LsZI"
      },
      "source": [
        "<font color='firebrick'><b>Task D2: Investigate the steady states.</b></font><br>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "4EHN47qnLsZI"
      },
      "source": [
        "If all went well you obtained two distinct steady state values of P. This means the steady state the system eventually reaches depends on the initial concentration of P. To determine what initial conditions are required to end up in a particular steady state, we are going to plot and analyse the rate characteristics (rates on the y-axis against P on the x-axis). <br> In order to obtain the rate characteristics we will define functions, using python. These functions should return the rate (r_form or r_deg) for any given concentration of P. Below such a function is defined for $r_{form}$:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "ZtJklUIaLsZJ",
        "colab": {}
      },
      "source": [
        "def calc_r_form(P): # We provide the function with a value for P as an input, see examples below\n",
        "    r_form = 3 + 20/(1+ np.float_power(2,(-3*(P-5))))\n",
        "    \n",
        "    return r_form"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "-CvOke0jLsZK"
      },
      "source": [
        "We can now call the function calc_r_form and use it to calculate r_form for a given P, for example:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "7XTVqzmLLsZK",
        "colab": {}
      },
      "source": [
        "calc_r_form(0) # calculate r_form when P = 0 mM"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "c_1MbKMZLsZO",
        "colab": {}
      },
      "source": [
        "calc_r_form(6) # calculate r_form when P = 6 mM"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "ErkLOnJzLsZP"
      },
      "source": [
        "Below a similar function needs to be defined for r_deg. <font color='firebrick'> **Task D2.1 Please complete the function below:** [1 points]</b></font>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "O6g_7oD1LsZP",
        "colab": {}
      },
      "source": [
        "def calc_r_deg(P)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "nWpFxk2wLsZR"
      },
      "source": [
        "The two functions above can now be used to obtain the rates between P = 0 mM and P = 11 mM (or any other range). It would be tedious to write a new line of code for each different value of P we would like the rates for:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "iu58koo92PCi",
        "colab": {}
      },
      "source": [
        "calc_r_deg(0)\n",
        "calc_r_deg(0.1)\n",
        "calc_r_deg(0.2)\n",
        "calc_r_deg(0.3)\n",
        "# etc...."
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "zz9s-mlE2Svx"
      },
      "source": [
        "Instead of writing many lines, we will use a so-called numpy array. Such an array holds multiple values. Below such an array containing P concentrations between 0 and 11 mM is defined."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "g6HsY4cXLsZR",
        "colab": {}
      },
      "source": [
        "P_concentrations = np.linspace(0,11,400)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "JRugiDieLsZS",
        "colab": {}
      },
      "source": [
        "P_concentrations"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "knX0adD7LsZU"
      },
      "source": [
        "As can been seen above, P_concentrations now contains 400 different concentrations of P, all between 0 and 11 (with steps of 11/400 between the values). Individual concentrations can be accesed as follows:"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "pJ1zXbK1LsZU",
        "colab": {}
      },
      "source": [
        "P_concentrations[0] # The first concentration stored in P_concentrations (Python starts counting at 0, not at 1!)"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "7PW-5S_tLsZW",
        "colab": {}
      },
      "source": [
        "P_concentrations[1]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "zyrGs2wVLsZZ",
        "colab": {}
      },
      "source": [
        "P_concentrations[399]"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "hR8S-fgnLsZb"
      },
      "source": [
        "<font color='firebrick'>**Task D2.2: Use the two functions (calc_r_form and calc_r_deg) to obtain two new numpy arrays. These new arrays should respectively contain the values of r_form and r_deg for the concentrations of P in P_concentrations. For example, r_form_array[0] should contain the value of r_form corresponding to the P concentration in P_concentrations[0]    [1 points] </b></font>. <br> Note that the functions can be applied to an array instead of a single concentration of P. You therefore do not need to use the function np.linspace to create these new arrays.**   \n"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "0YrN4sbfLsZb",
        "colab": {}
      },
      "source": [
        "# create the numpy array for r_form\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "LL3AFItPLsZd",
        "colab": {}
      },
      "source": [
        "# create the numpy array for r_deg\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "VqAM76nlLsZg"
      },
      "source": [
        "<font color='firebrick'> **Task D2.3: Next, we can use the three arrays (containing the P concentrations, rates for $r_{form}$, rates for $r_{deg}$) to plot the rate characteristics (y-axis: both rates in $mM \\; min^{-1}$, x-axis: [P] in mM). Please do so below using matplotlib (introduced above). [2 point]**</b></font>"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "T4qqRbWTLsZh",
        "colab": {}
      },
      "source": [
        "# PLOT FROM P = 0 mM to P = 11 mM"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "1yg5kdF_LsZj"
      },
      "source": [
        "<font color='firebrick'>Questions task D2  </b></font>: <br>\n",
        "    1. How can steady states be identified in the rate characteristic plotted above? [3 points]\n",
        "    2. Give the possible steady states of the system [1 point]\n",
        "    3. Under what initial P concentrations (if any) is each steady state reached? Use the rate characteristic plot to determine what the requirements are to reach each steady state. Check whether you are correct by running additional simulations below for different initial concentrations (you may add additional cells to the notebook below). [6 points]\n",
        "    \n",
        "   "
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "UZT3cD5oLaS0"
      },
      "source": [
        "Answers:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "2U5Bwm990I5r"
      },
      "source": [
        "<h2 align=\"left\"><font color='navy'> E. Analyse the model of blood coagulation pathway[11 points]</font></h2> \n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "llZyC2HJ0I5s"
      },
      "source": [
        "Up untill now we wrote down all models in this tutorial ourselves. In this exercise we will import a fully functional model instead. As mentioned in the lecture, there are databases with existing published models. It's imporant to write clear models to increase reproducibility of research.\n",
        "In this exercies we will import the model of blood coagulation pathway from BioModels database and perform some analysis.\n",
        "\n",
        "Coagulation is the process by which a blood clot is formed in order to stop bleeding. In a series of reactions, thrombin is activated. Thrombin is the the most important constituent of the coagulation cascade in terms of its activation role in the blood clotting process.\n",
        "Usually, the coagulation process is under the inhibitory control of several inhibitors that limit the clot formation (here, TFPI and ATIII)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "cy9vLVHr0I5t"
      },
      "source": [
        "<img src=\"https://github.com/A2P2/Tools-of-Sys-Bio_RUG/blob/master/Pictures/Hockin2002_fig.png?raw=true\" align=\"center\" width=\"250\">\n",
        "\n",
        "<center> <b>Figure 3:</b> Simplistic representation of blood coagulation pathway </center>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "LK-VPaZx0I5t"
      },
      "source": [
        "<font color='firebrick'><b>Task E1: Export the model from BioModels [2 points].</b></font><br>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "gB-6572z0I5t"
      },
      "source": [
        "To find the model go to https://www.ebi.ac.uk/biomodels/. There look in the search line for 'Hockin2002_BloodCoagulation'. In the Files tab find and copy the download link to the .xml file of the model (Right mouse button -> Copy Link Location, exact wording depends on the browser that you use). Copy paste the link in the line of code below."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "E6UPdPSu0I5t",
        "colab": {}
      },
      "source": [
        "model = te.loadSBMLModel(\"PasteTheLinkHereBetweenQuotes\")  "
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "5JX7-oKj0I5v",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE MODEL HERE\n",
        "#SIMULATE FROM START=0 TO END=1000\n",
        "#IN THIS MODEL UNIT FOR TIME IS SECONDS AND FOR CONCENTRATION IS M\n",
        "\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "Llxj05XG0I5w"
      },
      "source": [
        "<font color='firebrick'><b>Task E2: Plot only time trace of thrombin (IIa) [2 points].</b></font><br>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "t6qA-Ppr0I5x"
      },
      "source": [
        "As you see in your graph, the model has a lot of components. To simplify the graph, let's simulate only behaviour of thrombin (called IIa in the model). To do so you need to pass 'selections' to Tellurium:``` model.simulate(start=0, end=1000, points=2000, selections=['time','IIa'])``` And then continue with plotting."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "Pd-vCEw60I5x",
        "colab": {}
      },
      "source": [
        "# SIMULATE AND PLOT HERE"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "VdLjr7AS0I5z"
      },
      "source": [
        "Let's now see how the model string looks like. Execute the next line and scroll the output."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "ekLiWTjt0I5z",
        "colab": {}
      },
      "source": [
        "model.getAntimony().split('\\n')"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "03OFNojj0I52"
      },
      "source": [
        "<font color='firebrick'><i> Question E2.1: \tWhat are the default initial concentrations of TFPI and ATIII in the model? [1 points]</i></font>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "CvCnK22rogst"
      },
      "source": [
        "Answer E2.1:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "QrIzznaW0I52"
      },
      "source": [
        "<font color='firebrick'><b>Task E3: Explore influence of inhibitors on the cascade. [3 points]</b></font><br>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "jfMVpWiO0I52"
      },
      "source": [
        "For this task, separetely change the initial concentrations of both inhibitors to 0. Plot the results. Describe how the absence of each inhibitor affects the shape and maximal value of the thrombin concentration curve."
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "XozMZUsQ0I52",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE MODEL WITH TFPI=0 HERE\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "code",
      "metadata": {
        "colab_type": "code",
        "id": "_8xCn30-0I55",
        "colab": {}
      },
      "source": [
        "#SIMULATE AND PLOT THE MODEL WITH ATIII=0 HERE\n"
      ],
      "execution_count": null,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "srGYvWrx0I57"
      },
      "source": [
        "<font color='firebrick'><i> Question E3.1:  Why should the coagulation cascade be inhibited? What could the physiological consequence(s) if the coagulation cascade is not inhibited? [3 points]</i></font>\n"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "LuNk7I-U1GaS"
      },
      "source": [
        "Answer E3.1:"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "W9cNb5nJch13"
      },
      "source": [
        "### References ###"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "colab_type": "text",
        "id": "suqwiUJwch14"
      },
      "source": [
        "1. Practicum *Kinetic models* (2016) T. Kimkes, V. Takhaveev, M. Heinemann\n",
        "2. Tellurium web-site: https://tellurium.readthedocs.io/en/latest/installation.html\n",
        "3. Porubsky, V. Tellurium tutorial. https://github.com/vporubsky/tellurium-libroadrunner-tutorial\n",
        "4. Ferrell, J. E., Tsai, T. Y. C., & Yang, Q. (2011). Modeling the cell cycle: why do certain circuits oscillate?. Cell, 144(6), 874-885.\n",
        "5. Fung, E., Wong, W. W., Suen, J. K., Bulter, T., Lee, S. G., & Liao, J. C. (2005). A synthetic gene–metabolic oscillator. Nature, 435(7038), 118-122.\n",
        "6. Hockin MF, Jones KC, Everse SJ, Mann KG. A model for the stoichiometric regulation of blood coagulation. J Biol Chem. 2002 May 24;277(21):18322-33."
      ]
    }
  ]
}