2-ode-simulation.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Simulation of ODEs in Julia"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To simulate (systems of) ODEs in Julia, we use the [DifferentialEquations](https://docs.sciml.ai/DiffEqDocs/stable/) module. This makes it easy to define all kinds of ODE systems, from very simple towards increasingly complicated definitions. This overview shows some of the basic functionality of ODE systems, from simulating a basic set of equations to the use of callbacks to perturb the system at given timepoints. We start with importing some of the required modules into Julia with the `using` keyword. In this example, we only use the `DifferentialEquations` module and the `Plots` module for plotting our final solutions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "using DifferentialEquations, Plots, Sundials"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Lotka-Volterra System\n",
    "A simple example of an ODE system that can easily be simulated is are the Lotka-Volterra equations. These are a classic example of a predator-prey system, or a feast-famine regime. The equations are as follows.\n",
    "$$\n",
    "\\begin{align}\n",
    "\\dot{x} &= \\alpha x - \\beta x y \\\\\n",
    "\\dot{y} &= \\gamma x y - \\delta y\n",
    "\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "The system contains four parameters, depicted by $[\\alpha, \\beta, \\gamma, \\delta]$. Defining a system of ODEs in Julia is similar to matlab. However, even though it is not necessary for the solver to work, we add a parameter `du` to the input of our ODE, which makes the function a *mutating* function (indicated by the exclamation mark at the end of the function name, which is a Julia convention). Instead of defining a new vector of derivatives each function call, we mutate the input derivative vector, which is beneficial for speed and memory usage.\n",
    "\n",
    "Following the system definition, we define an initial condition and a timespan for our solver. We finally give values to our parameter vector, before we can ask Julia to define an `ODEProblem`. This is an object defined by the `DifferentialEquations.jl` library that works with its `solve` function, and optionally performs some optimizations under the hood. Even though `solve` also has function arguments, we leave it to its simplest form here by just calling it on our recently defined `ODEProblem`.\n",
    "\n",
    "Plotting our solution is then just a small step away. We have defined some parameters in the plot function to make it look nicer, but `plot(solution)` would have worked as well."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/Users/max/TUe/lotka-volterra.svg\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "function lotka_volterra!(du, u, parameters, t)\n",
    "    du[1] = parameters[1].*u[1] .- parameters[2].*u[1].*u[2]\n",
    "    du[2] = parameters[3].*u[1].*u[2] .- parameters[4].*u[2]\n",
    "    nothing\n",
    "end\n",
    "\n",
    "# initial condition and time span\n",
    "u0 = [10., 1.]\n",
    "tspan = (0., 80.)\n",
    "\n",
    "# parameters\n",
    "parameters = [0.8, 0.4, 0.1, 0.3]\n",
    "\n",
    "# simulation\n",
    "problem = ODEProblem(lotka_volterra!, u0, tspan, parameters)\n",
    "solution = solve(problem)\n",
    "\n",
    "plot(solution, labels = [\"Prey\" \"Predator\"], linewidth=2, title=\"Example Lotka-Volterra System\", xlabel=\"Time\", ylabel=\"Population Size\")\n",
    "savefig(\"../lotka-volterra.svg\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chemical Reaction System\n",
    "\n",
    "A different system could be a simple homeostatic system, influenced by possible perturbations $X$. We can define the system as follows.\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "\\dot{X} &= -p_1 X \\\\\n",
    "\\dot{A} &= p_1 X + \\frac{p_2}{1+C} - k_1A - c_1A \\\\\n",
    "\\dot{B} &= k_1A - k_2B \\\\\n",
    "\\dot{C} &= \\frac{k_2B}{1+C} - c_2C  \\\\\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "Implementing this system in Julia is highly similar to the Lotka-Volterra system. We first set the initial perturbation to zero, to explore the steady state evolution of the system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/Users/max/TUe/steady-state.svg\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "function chemsys!(du, u, p, t)\n",
    "    X, A, B, C = u\n",
    "    p1, p2, k1, k2, c1, c2 = p\n",
    "    du[1] = -p1.*X\n",
    "    du[2] = p1.*X+((p2)./(1 .+C)) .- (k1+c1).*A\n",
    "    du[3] = k1.*A .- k2.*B\n",
    "    du[4] = (k2.*B)./(1 .+C) .- c2.*C\n",
    "end;\n",
    "\n",
    "# parameters\n",
    "u0 = [0.2, 1., 0., 0.]\n",
    "p = [0.1, 0.1, 0.15, 0.07, 0.01, 0.2]\n",
    "tspan = (0., 150.)\n",
    "\n",
    "# and create an ode problem to solve\n",
    "problem = ODEProblem(chemsys!, u0, tspan, p);\n",
    "\n",
    "# the solution can be obtained easily\n",
    "solution = solve(problem);\n",
    "\n",
    "plot(solution, linewidth=2, labels=[\"Input\" \"A\" \"B\" \"C\"], xlabel=\"Time\", ylabel=\"Concentration\", title=\"Steady-State Evolution of Chemical System\")\n",
    "savefig(\"../steady-state.svg\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As can be seen, the system moves to a steady state quite quickly. We can extend the timespan and add some perturbations. In this case, we have decided to add three perturbations at $t = 100, t = 200, t = 300$, with slightly increased peak heights. We can model this as a callback function, and add the times of the dose to the `tstops` argument to make sure that the integrator stops at the times of the perturbations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\"/Users/max/TUe/julia-ode-tutorial/perturbed-system.svg\""
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# change the timespan\n",
    "problem = remake(problem; tspan=(0.,500.))\n",
    "\n",
    "# set the times of the dose\n",
    "dosetimes = [100. 200. 300.]\n",
    "\n",
    "# set the callback condition (if t is in the dosetimes array, apply a perturbation)\n",
    "condition(u,t,integrator) = t ∈ dosetimes\n",
    "\n",
    "# function to determine the perturbation size\n",
    "affect! = function(integrator) \n",
    "    integrator.u[1] += 0.4.*integrator.t/100\n",
    "end\n",
    "\n",
    "# combine the condition and the perturbation into a callback\n",
    "cb = DiscreteCallback(condition,affect!)\n",
    "\n",
    "# compute the solution with perturbation\n",
    "solution_with_perturbation = solve(problem,Tsit5(),callback=cb, tstops=dosetimes)\n",
    "\n",
    "# plot the result\n",
    "plot(solution_with_perturbation, labels= [\"Input\" \"A\" \"B\" \"C\"], linewidth=2, title=\"Perturbed System\", xlabel=\"Time\", ylabel=\"Concentration\")\n",
    "\n",
    "savefig(\"perturbed-system.svg\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can see that the function restores to the steady-state solution after the perturbations. "
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.8.2",
   "language": "julia",
   "name": "julia-1.8"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.8.2"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "b0fa6594d8f4cbf19f97940f81e996739fb7646882a419484c72d19e05852a7e"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}