{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Another Delayed Nonlinear Negative Feedback Oscillator\n",
    "\n",
    "Implement the 7-state model of Becker-Weimann (Becker-Weimann, Sabine, Jana Wolf, Hanspeter Herzel, and Achim Kramer. 2004. “Modeling Feedback Loops of the Mammalian Circadian Oscillator.” Biophysical Journal 87 (5): 3023–34. https://doi.org/10.1529/biophysj.104.040824.).\n",
    "\n",
    "The parameters are\n",
    "- v1b : maximal rate of Per2/Cry transcription\n",
    "- k1b : activation threshold for Per2/Cry transcription\n",
    "- k1i : inhibition constant of Per2/Cry transcription\n",
    "- c: concentration of constituative activator\n",
    "- p: Hill coefficient of inhibition of Per2/Cry transcription\n",
    "- k1d: degradation rate of Per2/Cry mRNA\n",
    "- k2b: complex formation rate of Per2/Cry mRNA\n",
    "- q: Number of PER2/CRY complex forming subunits\n",
    "- k2d: degradation rate of cytoplasmic PER2/CRY\n",
    "- k2t: nuclear import rate of teh PER2/CRY complex\n",
    "- k3t: nuclear export rate of teh PER2/CRY complex\n",
    "- k3d: degradation rate of nuclear PER2/CRY\n",
    "- v4b: maximal rate of Bmal1 transcription\n",
    "- k4b : activation threshold for Bmal1 transcription\n",
    "- r: Hill coefficient of activation of Bmal1 transcription\n",
    "- k4d: degradation rate of Bmal1 mRNA\n",
    "- k5b: translation rate of BMAL1\n",
    "- k5d: degradation rate of cytoplasmic BMAL1\n",
    "- k5t: nuclear import rate of BMAL1\n",
    "- k6t: nuclear export rate of BMAL1\n",
    "- k6d: degradation rate of nuclear BMAL1\n",
    "- k6a: activation rate of nuclear BMAL1\n",
    "- k7a: deactivation rate of nuclear BMAL1\n",
    "- k7d: degradataion rate of nuclear BMAL1\n",
    "\n",
    "The model is\n",
    "\\begin{eqnarray*}\n",
    "\\frac{dy1}{dt} &=& f_1 - k_{1d}\\cdot y1\\\\\n",
    "\\frac{dy2}{dt} &=& k_{2b}\\cdot y1^q - k_{2d}\\cdot y2 - k_{2t}\\cdot y2 + k_{3t}\\cdot y3\\\\\n",
    "\\frac{dy3}{dt} &=& k_{2t}\\cdot y2 - k_{3t}\\cdot y3 - k_{3d}\\cdot y3 \\\\\n",
    "\\frac{dy4}{dt} &=& f_2 - k_{4d}\\cdot y4\\\\\n",
    "\\frac{dy5}{dt} &=& k_{5b}\\cdot y4 - k_{5d}\\cdot y5 - k_{5t}\\cdot y5 + k_{6t}\\cdot y6 \\\\\n",
    "\\frac{dy6}{dt} &=& k_{5t}\\cdot y5 - k_{6t}\\cdot y6 - k_{6d}\\cdot y6 + k_{7a}\\cdot y7 - k_{6a}\\cdot y6\\\\\n",
    "\\frac{dy7}{dt} &=& k_{6a}\\cdot y6 - k_{7a}\\cdot y7 - k_{7d}\\cdot y7\\\\\n",
    "\\end{eqnarray*}\n",
    "\n",
    "where\n",
    "\n",
    "\\begin{equation*}\n",
    "f_1 = \\frac{v_{1b}\\cdot(y7+c)}{k_{1b}\\cdot(1+(y3/k_{1i})^p)+(y7+c)}\n",
    "\\end{equation*}\n",
    " \n",
    "and\n",
    "\n",
    "\\begin{equation*}\n",
    "f_2 = \\frac{v_{4b}\\cdot y3^r}{k_{4b}^r+y3^r}\n",
    "\\end{equation*}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "def bw_oscillator( t, statevars, params ):\n",
    "    y1 = statevars[0];\n",
    "    y2 = statevars[1];\n",
    "    y3 = statevars[2];\n",
    "    y4 = statevars[3];\n",
    "    y5 = statevars[4];\n",
    "    y6 = statevars[5];\n",
    "    y7 = statevars[6];\n",
    "\n",
    "    v1b = params[0] # maximal rate of Per2/Cry transcription\n",
    "    k1b = params[1] # activation threshold for Per2/Cry transcription\n",
    "    k1i = params[2] # inhibition constant of Per2/Cry transcription\n",
    "    c = params[3] # concentration of constituative activator\n",
    "    p = params[4] # Hill coefficient of inhibition of Per2/Cry transcription\n",
    "    k1d = params[5] # degradation rate of Per2/Cry mRNA\n",
    "    k2b = params[6] # complex formation rate of Per2/Cry mRNA\n",
    "    q = params[7] # Number of PER2/CRY complex forming subunits\n",
    "    k2d = params[8] # degradation rate of cytoplasmic PER2/CRY\n",
    "    k2t = params[9] # nuclear import rate of teh PER2/CRY complex\n",
    "    k3t = params[10] # nuclear export rate of teh PER2/CRY complex\n",
    "    k3d = params[11] # degradation rate of nuclear PER2/CRY\n",
    "    v4b = params[12] # maximal rate of Bmal1 transcription\n",
    "    k4b = params[13] # activation threshold for Bmal1 transcription\n",
    "    r = params[14] # Hill coefficient of activation of Bmal1 transcription\n",
    "    k4d = params[15] # degradation rate of Bmal1 mRNA\n",
    "    k5b = params[16] # translation rate of BMAL1\n",
    "    k5d = params[17] # degradation rate of cytoplasmic BMAL1\n",
    "    k5t = params[18] # nuclear import rate of BMAL1\n",
    "    k6t = params[19] # nuclear export rate of BMAL1\n",
    "    k6d = params[20] # degradation rate of nuclear BMAL1\n",
    "    k6a = params[21] # activation rate of nuclear BMAL1\n",
    "    k7a = params[22] # deactivation rate of nuclear BMAL1\n",
    "    k7d = params[23] # degradataion rate of nuclear BMAL1\n",
    "        \n",
    "    dydt = np.zeros( (7,) )\n",
    "    f1 = \n",
    "    f2 = \n",
    "\n",
    "    dydt[0] = \n",
    "    dydt[1] = \n",
    "    dydt[2] = \n",
    "    dydt[3] = \n",
    "    dydt[4] = \n",
    "    dydt[5] = \n",
    "    dydt[6] = \n",
    "\n",
    "    return dydt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Simulate the model\n",
    "\n",
    "Simulate the model with the parameters given in Table of the paper. I have already coded them up, because there are so many!\n",
    "\n",
    "We will assume the model reaches its limit cycle within 200 hours of simulation. Simulate the model first for 200 hours. Then start a new simulation from the end time point of the first simulation. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.integrate\n",
    "\n",
    "v1b = 9 # maximal rate of Per2/Cry transcription\n",
    "k1b = 1 # activation threshold for Per2/Cry transcription\n",
    "k1i = 0.56 # inhibition constant of Per2/Cry transcription\n",
    "c = 0.01 # concentration of constituative activator\n",
    "p = 8 # Hill coefficient of inhibition of Per2/Cry transcription\n",
    "k1d = 0.12 # degradation rate of Per2/Cry mRNA\n",
    "k2b = 0.3 # complex formation rate of Per2/Cry mRNA\n",
    "q = 2 # Number of PER2/CRY complex forming subunits\n",
    "k2d = 0.05 # degradation rate of cytoplasmic PER2/CRY\n",
    "k2t = 0.24 # nuclear import rate of teh PER2/CRY complex\n",
    "k3t = 0.02 # nuclear export rate of teh PER2/CRY complex\n",
    "k3d = 0.12 # degradation rate of nuclear PER2/CRY\n",
    "v4b = 3.6 # maximal rate of Bmal1 transcription\n",
    "k4b = 2.16 # activation threshold for Bmal1 transcription\n",
    "r = 3 # Hill coefficient of activation of Bmal1 transcription\n",
    "k4d = 0.75 # degradation rate of Bmal1 mRNA\n",
    "k5b = 0.24 # translation rate of BMAL1\n",
    "k5d = 0.06 # degradation rate of cytoplasmic BMAL1\n",
    "k5t = 0.45 # nuclear import rate of BMAL1\n",
    "k6t = 0.06 # nuclear export rate of BMAL1\n",
    "k6d = 0.12 # degradation rate of nuclear BMAL1\n",
    "k6a = 0.09 # activation rate of nuclear BMAL1\n",
    "k7a = 0.003 # deactivation rate of nuclear BMAL1\n",
    "k7d = 0.09 # degradataion rate of nuclear BMAL1\n",
    "\n",
    "params = (v1b,k1b,k1i,c,p,k1d,k2b,q,k2d,k2t,k3t,k3d,v4b,k4b,r,k4d,k5b,k5d,k5t,k6t,k6d,k6a,k7a,k7d)\n",
    "\n",
    "mf = lambda t,y : bw_oscillator(t,y, params)\n",
    "\n",
    "sol1 = scipy.integrate.solve_ivp( fun=mf, t_span=(0,200), y0=0.1*np.ones((7,)), \n",
    "                                t_eval=np.arange(0,200,.1) )\n",
    "sol2 = scipy.integrate.solve_ivp( fun=mf, t_span=(0,100), y0=sol1.y[:,-1], \n",
    "                                t_eval=np.arange(0,100,.1) )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot the simulation\n",
    "Plot both simulations in separate subplots. The first should show the transient (the path it takes to reach the oscillation's limit cycle). The second should show the oscillation on its limit cycle)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Get to know the regulation function for Per2/Cry\n",
    "How does Per2/Cry regulation (f1) change over time and how does it relate to per2/cry mRNA production? Compute f1 using the sol2. Then plot Per2/Cry mRNA over time (sol2.y[0,:]), and f 1 over time in separate subplots."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# write code here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "How the activation function depend on to y7 and y3? Plot a heatmap of the activation function over the range of values (N=100) of y3 and y7 from the simulation. Then, plot the limit cycle on top.\n",
    "\n",
    "The activation for each value of y3 and y7 will be stored in a 100x100 matrix. To plot the heatmap, we will need to use plt.imshow, which is meant to show images. The y-axis is the rows and is numbered from the top to the bottom. The x-axis is the columns and is numbered from left to right. To plot the limit cycle, we will need the indices into the y3 and y7 ranges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make arrays of 100 evenly spaced values covering the ranges of y7 and y4\n",
    "\n",
    "# Make the matrix. 1 row for each value of y3, 1 column for each value of y7\n",
    "        \n",
    "# plot the heatmap.\n",
    "# Since y7 changes across columns, y7 indices are shown along the x-axis\n",
    "# Since y3 changes across rows, y3 indices are shown along the y-axis\n",
    "\n",
    "# At each time step, determine which value in our list of 100 y7 values\n",
    "# is closest to that value at that time step. Use that index to determine where\n",
    "# to plot the point. Do the same for y3.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Find CT0\n",
    "\n",
    "Circadian Time (CT) is a standardardized system for indicating what \"time\" the clock is telling. Each day goes from CT0 to CT24 (well, sort, of CT0=CT24). CT0 refers to the internal time when we expect dawn. CT12 refers to the internal time when we expect dusk. It has been observed that Per2 mRNA peaks sometime between CT6 and CT10. For our model, we will find the peak of Per2 mRNA, and say CT0 is 7 hours before that peak. (Note that circadian hours and hours are not exactly equal because 1 circadian hour is 1/24th of a cycle. If the period of the cycle isn't exactly 24 h, then they are not equivalent. But given that our goal is simply to get that peak in the range of CT6 to CT10, this is good enough.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look for the peak of Per2/Cry mRNA in the last 24h of the simulation.\n",
    "\n",
    "# Now, run the simulation from CT0 for 100 hours.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}