{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "msxy7km-lwO9"
   },
   "source": [
    "### **Course**: BIO-341 [_Dynamical systems in biology_](https://moodle.epfl.ch/course/info.php?id=14291)\n",
    "\n",
    "**Professor**: _Julian Shillcock_ & _Felix Naef_\n",
    "\n",
    "SSV, BA5, 2025"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "03GfzrTflwPA"
   },
   "outputs": [],
   "source": [
    "#import important libraries\n",
    "import numpy as np\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from IPython.display import set_matplotlib_formats\n",
    "from scipy.integrate import odeint\n",
    "import math as mt\n",
    "from mpl_toolkits.mplot3d import Axes3D\n",
    "from matplotlib.lines import Line2D\n",
    "import random"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Circadian oscillators"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Sample exam question: Coupled oscillators (Paper and pencil)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The phases of these two oscillators with natural frequency (fréquence propre) $\\omega_1$ and $\\omega_2$ are written $\\phi_1(t)$ and\n",
    "$\\phi_2(t)$. \n",
    "\n",
    "They are coupled according to the following model:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\frac{d\\phi_1}{dt} = \\omega_1 + K f(\\phi_2 - \\phi_1)\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{d\\phi_2}{dt} = \\omega_2  + K f(\\phi_1 - \\phi_2)\n",
    "\\end{equation}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with $f(x) = \\sin(2x)$ "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Determine whether the two phases synchronise, i.e. whether the phase difference $\\alpha =\\phi_2-\\phi_1$ reaches a stable fixed point."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Form the differential equation for the phase difference $\\alpha =\\phi_2-\\phi_1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Explain whether this model allows the two phases to be synchronised, and if so, under what conditions.\n",
    "conditions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. What happens to $\\alpha$ when $K$ goes to infinity?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "K3ZLuiG6lwPB"
   },
   "source": [
    "## 2. Circadian oscillators of a generic non-linear oscillator (Python)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "x20J1jcklwPC"
   },
   "source": [
    "Here, we first simulate a toy model of the circadian clock consisting of three variables that implements the now famous negative feedback loop oscillator (See the 2017 Nobel price in Physiology and Medicine)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Cgy3UKkQlwPC"
   },
   "source": [
    "## A three-variable  model of a circadian oscillator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7_A7p_uylwPC"
   },
   "source": [
    "This is a highly simplified model of a circadian oscillator (see Feedback of the Drosophila period gene product on circadian cycling of its messenger RNA levels, Hardin P, Hall JC and Rosbash M, Nature 1990). It takes into account some basic ingredients, notably the negative feedback loop. A clock gene mRNA (X) produces a clock protein (Y) which, in turn, activates a transcriptional inhibitor (Z) of gene X."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3xLKb7X-lwPD"
   },
   "source": [
    "\\begin{equation}\n",
    "\\frac{dX}{dt} = v_{1} \\frac{K_{1}^{4}}{K_{1}^{4} + Z^ {4}} - v_{2} \\frac{X}{K_{2} + X} \n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dY}{dt} = k_{3}X  - v_{4} \\frac{Y}{K_{4} + Y} \n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dZ}{dt} = k_{5}Y  - v_{6} \\frac{Z}{K_{6} + Z}  \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Z5e3DodrlwPD"
   },
   "source": [
    "### The Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "7pyvpD5rlwPE"
   },
   "source": [
    "1) Explain the different terms in the equations, and all the parameters. In particular explain the terms containing the fractions."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "s9xuhL3BlwPE"
   },
   "source": [
    "2) Based on the literature, discuss plausible genes/proteins that could represent the $X,Y,Z$ variables."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "1VJ1_eV-lwPF"
   },
   "source": [
    "3) Discuss/criticize the main assumptions of the model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zvj97huTlwPF"
   },
   "source": [
    "### Simulation of the Model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "X42MRGVTlwPG"
   },
   "source": [
    "4) Using the values $v_{1} = 0.7 \\,nMh^{-1}, v_{2} = v_{4} = v_{6} = 0.35\\, nMh^{-1}, K_{1}=K_{2}=K_{4}=K_{6} = 1\\,nM$ and $k_{3}=k_{5}=0.7\\,h^{-1}$, simulate the model:  set appropriate initial conditions and time integration parameters to obtain a limit cycle. Plot some representative X(t), Y(t), and/or Z(t), using appropriate initial conditions. Finally, plot the trajectory in the X, Y and Z phase portrait (using `ax = plt.axes(projection='3d')`).\n",
    "   \n",
    ">Hint: to visualise the trajectories, plot only the last 5 % of the solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "-dMtmoyUlwPI"
   },
   "outputs": [],
   "source": [
    "#Solving the differential equations\n",
    "\n",
    "Tmax=5000\n",
    "dt=0.01\n",
    "tspan = np.arange(0, Tmax, dt)\n",
    "\n",
    "# Initial conditions setting \n",
    "X0=[0.14,0.18,1.8]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Definition of the model \n",
    "\n",
    "def model(s, t):\n",
    "    \n",
    "    # add here the equations of the model\n",
    "    \n",
    "    return x_dot, y_dot, z_dot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "8CgzhT1jlwPJ"
   },
   "source": [
    "5) Use Period_finder on a long (many periods) trajectory on x. \n",
    "Comment the period distribution you find and its point estimate.\n",
    "What happens if you use the y or z trajectory to evaluate the period?\n",
    "\n",
    ">The Period_finder function takes as an imput the x (or y or z) vector of coordinates and a vector of equally spaced times at wich the given coordinate was obtained with a simulation. \n",
    "It takes the x limits of the plot as a possible input, with default values of 0 to 50.\n",
    "It returns a plot of the period distribution and a point estimate for the period.\n",
    "If you want to understand how this function was built, read \"A. few words on the discrete Fourier transform\".\n",
    "\n",
    "> Hint: use the last 20% of the solution to find the period"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "27q3ai3MnYnH"
   },
   "outputs": [],
   "source": [
    "def Period_finder(x, tspan, xlim=[0,50], ToPrint=True):\n",
    "    \n",
    "    signal = x\n",
    "    \n",
    "    omega = np.fft.fft(signal)\n",
    "    \n",
    "    modes = np.arange(omega.size)\n",
    "    \n",
    "    t_dist = (tspan[-1]-tspan[0])\n",
    "    \n",
    "    omega_cut = len(omega)//2\n",
    "    \n",
    "    modes = modes[1:omega_cut]\n",
    "    omega = omega[1:omega_cut]\n",
    "    \n",
    "    periods = t_dist/modes\n",
    "    \n",
    "    abs_o = np.absolute(omega)\n",
    "    max_o = np.argmax(abs_o)\n",
    "    \n",
    "    period_estimate = periods[max_o]\n",
    "    \n",
    "    if ToPrint:\n",
    "        plot_period(period_estimate, periods, abs_o)\n",
    "        print(\"the period value is about \" + str(round(period_estimate,2)) + \" hours\")\n",
    "    \n",
    "    return (period_estimate, periods, abs_o)\n",
    "\n",
    "def plot_period(period_estimate, periods, abs_o,xlim=[0,50]):\n",
    "    fig,ax= plt.subplots()\n",
    "\n",
    "    ax.axvline(period_estimate, ls='--', c='k')\n",
    "    \n",
    "    ax.plot(periods,abs_o/np.sum(abs_o))\n",
    "    \n",
    "    ax.set_xlim(xlim)\n",
    "    ax.set_xlabel(\"period\")\n",
    "    ax.set_ylabel(\"propobabilty density\")\n",
    "    \n",
    "    return (fig ,ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "3On_-PGWlwPM"
   },
   "source": [
    "#### Hopf bifurcation.\n",
    "6) Vary the value of the transcription rate $v_{1}$ in the interval $(0,5]$ $nMh^{-1}$. You can plot some representative trajectories (see the code before and replace the $v_{1}$). Plot and discuss the bifurcation diagram (show $X_{min}$ and $X_{max}$ in function of $v_{1}\\in\\left(0,5\\right]$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sample_frac=0.05\n",
    "\n",
    "def model(s, t, v):\n",
    "    # define the model\n",
    "    return x_dot, y_dot, z_dot\n",
    "    \n",
    "f, axs = plt.subplots(4,5, figsize=(20,10))\n",
    "axs=axs.flatten()\n",
    "cmap = matplotlib.cm.viridis\n",
    "\n",
    "# We define a range in which we vary v\n",
    "vspan=np.linspace(0.1,5,20)\n",
    "\n",
    "lastpart=int(len(tspan)*samp_frac) \n",
    "\n",
    "X_lims=[]\n",
    "for n,v in enumerate(vspan):\n",
    "    # solve the system\n",
    "    sol=odeint(model, X0, tspan, args=(v,))\n",
    "    # plot representative trajectories x vs y        \n",
    "    \n",
    "X_lims=np.array(X_lims)\n",
    "\n",
    "# Bifurcation diagram\n",
    "# plot the minimum and the maximum of the solution according to v"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Entrainement of a generic non-linear oscillator (Python)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take the model of exercise 2 and add an external periodic entrainement with period $T$ as follows:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "\\frac{dX}{dt} =  v_{1\\_ent}(t) * \\frac{K_{1}^{4}}{K_{1}^{4} + Z^ {4}} - v_{2} \\frac{X}{K_{2} + X}\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dY}{dt} = k_{3}X  - v_{4} \\frac{Y}{K_{4} + Y} \n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dZ}{dt} = k_{5}Y  - v_{6} \\frac{Z}{K_{6} + Z}  \n",
    "\\end{equation}\n",
    "\n",
    "where $v_{1\\_ent}(t) = v_{1} * A (1 + \\sin(\\frac{2\\pi}{T} t))$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Implement the entrainment in the model with $A = 0.01$ and $A = 0.2$. Is the oscilaltor stably entrained? If, what period do you expect then? Verify your prediction by simulation.\n",
    "   \n",
    ">Hint: to visualise the trajectories, plot only the last 5 % of the solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_external = 25\n",
    "Tmax=T_external * 100\n",
    "dt=0.01\n",
    "tspan = np.arange(0, Tmax, dt)\n",
    "# Initial conditions setting \n",
    "X0=[0.14,0.18,1.8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Definition of the model with entrainment\n",
    "def model(s, t, A, T_external=25):\n",
    "    # define the model here\n",
    "    return x_dot, y_dot, z_dot"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. Vary the value of the entrainment strength $A$ in the interval $[0,0.1]$. Which approximative value of A induces a change in period in the model?  You can plot representative trajectories and period diagram using the $period\\_finder$ provided function. Discuss your findings.\n",
    "\n",
    ">Hint: to calculate the period, use the last 20 % of the solution and to visualise the trajectories, plot only the last 5 % of the solution\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_external=25\n",
    "Tmax = T_external * 100 \n",
    "dt=0.01\n",
    "tspan = np.arange(0, Tmax, dt)\n",
    "\n",
    "# this tspan is used to calculate the period\n",
    "tspan_80percent = int(0.8*len(tspan))\n",
    "\n",
    "# this tspan is used for the plotting of the solution\n",
    "tspan_5percent = int(0.05*len(tspan))\n",
    "\n",
    "# Initial Conditions\n",
    "X0=[0.15, 0.3, 1.8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(s, t, A):\n",
    "# define your model here\n",
    "    return x_dot, y_dot, z_dot\n",
    "# Define the range in which A vary\n",
    "# A\n",
    "for n,a in enumerate(vspan):\n",
    "    # solve the system\n",
    "    sol=odeint(model, X0, tspan, args=(a,))\n",
    "# plot representative trajectories x vs y\n",
    "# plot the period estimate using the period_finder function"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. To analyse the change of period in relation to the entrainment strength $A$, simulate the model without entrainment for a time $T0$ and then switch on the entrainment. Vary the strength of the entrainment as in 3.2. Discuss the effect of the strength of the entrainment on the period (compare the period before and after switching the entrainment). Discuss you findings.\n",
    "\n",
    ">Hint: to calculate the period, here use the last 40 % of the solution and to visualise the trajectories, plot only the last 5 % of the solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "T_external=25\n",
    "\n",
    "# Total time of the simulation of the model\n",
    "Tmax = T_external * 100 #here we simulate a total of 100 cycles, you can also change this value\n",
    "dt=0.01\n",
    "\n",
    "# Time after the entrainment should start\n",
    "T0 = T_external * 40\n",
    "n_entrain= int(T0 / dt)\n",
    "\n",
    "# Initial Conditions\n",
    "X0=[0.15, 0.3, 1.8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def model(s, t, A):\n",
    "# define your model here\n",
    "    return x_dot, y_dot, z_dot\n",
    "# Define the range in which A vary\n",
    "# A\n",
    "# Define the two timespan vectors for before and after entrainment\n",
    "#\n",
    "for n,a in enumerate(vspan):\n",
    "# solve the system\n",
    "    sol=odeint(model, X0, tspan, args=(a,))\n",
    "# plot representative trajectories x vs y\n",
    "# plot the period estimate using the period_finder function"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "colab": {
   "collapsed_sections": [],
   "name": "Solutions 10.ipynb",
   "provenance": []
  },
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.11.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}