technical-notes/1-1-pacing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Technical note: Pacing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ODEs used in cardiac cellular electrophysiology are typically _forced_ by either applying a periodic stimulus current or by fixing the membrane potential to a stepwise constant value.\n",
    "In Myokit, this is called _pacing_.\n",
    "\n",
    "For efficient simulation, an _event-based pacing system_ is implemented in the Python classes `Protocol`, `ProtocolEvent`, and `PacingSystem`, and in the C header `pacing.h`.\n",
    "\n",
    "In addition, a _time-series pacing system_ is implemented for the CVODE single-cell simulation only.\n",
    "This is much less efficient, but may be useful in special cases.\n",
    "Time-series pacing is implemented in the C header `pacing.h`.\n",
    "\n",
    "For both systems, the output is a single valued, dimensionless quantity $x$, stored in the _pacing variable_.\n",
    "Finally, simulations can be run without pacing.\n",
    "Some notes on this are given in the section \"No pacing\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time\n",
    "\n",
    "The word _time_ in Myokit is used to denote a dimensionless clock variable  $t$, starting at $t=0$ and increasing indefinitely."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Event-based pacing\n",
    "\n",
    "Each event is defined as a 5-tuple $(x_i, t_i, d_i, p_i, m_i)$ or `(level, start, duration, period, multiplier)`.\n",
    "\n",
    "Here,\n",
    "- $x_i$ is a value for the pacing variable,\n",
    "- $t_i$ is the time the event starts and\n",
    "- $d_i$ is the event's duration.\n",
    "\n",
    "The values $p_i$ and $m_i$ are used to create periodic events.\n",
    "- When $p_i$ is set to any non-zero value, the event repeats every $p_i$ time units.\n",
    "- The event then occurs $m_i$ times, unless $m_i$ is set to zero in which case it recurs indefinitely."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Examples:\n",
    "\n",
    "1. A pulse of $x=1$ for $0.5$ time units, starting at $t=50$ and recurring indefinitely with a period of $1000$ time units is specified as `(1, 50, 0.5, 1000, 0)`.\n",
    "\n",
    "2. A singular pulse of $x=-40$, starting at $t=1000$ and lasting $100$ time units is specified as `(-40, 1000, 100, 0, 0)`.\n",
    "\n",
    "3. A periodic pulse of $p=1$ with duration $d_i=0.5$ that starts at $t=20$, $t=1020$ and $t=2020$ is specified as `(1, 20, 0.5, 1000, 3)`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Definitions\n",
    "\n",
    "- A singular event is an event with $p_i = 0$ and $m_i = 0$.\n",
    "- A singular event $i$ is active when $t_i \\leq t < t_i + d_i$.\n",
    "- A periodic event is an event with $p_i > 0$.\n",
    "- A periodic event with $m_i = 0$ occurs indefinitely.\n",
    "- An indefinitely occuring periodic event $i$ is active when $ t_i + k \\cdot r_i \\leq t < t_i + k \\cdot p_i + t_d$ where $k \\in \\{0,1,2,\\ldots\\}$.\n",
    "- A periodic event with $m_i > 0$ occurs a total of $m_i$ times.\n",
    "- A periodic event $i$ with $m_i > 0$ is active when $ t_i + k \\cdot r_i \\leq t < t_i + k \\cdot p_i + t_d$ where $k \\in \\{0,1,2,\\ldots,m_i-1\\}$.\n",
    "\n",
    "Restrictions:\n",
    "\n",
    "- All events with $p_i = 0$ must have $m_i = 0$.\n",
    "- The duration of a periodic event may never exceed its period: $d_i \\leq p_i$.\n",
    "\n",
    "<img src=\"figures/event.svg\" width=\"500px\" />"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Event-based pacing systems\n",
    "\n",
    "A system is an ordered sequence of events.\n",
    "It defines a pacing value $x$ and a special event, called the _current event_ that determines $x$'s value.\n",
    "If there is no current event, then $x=0$.\n",
    "\n",
    "Rules:\n",
    "\n",
    "- No two events in a pacing system may start or re-occur at the same time.\n",
    "- Whenever an event $i$ starts, it becomes the current event and the pacing variable's value changes to $x=x_i$.\n",
    "- When a second event $j$ starts during, or just after the current event, it becomes the current event and $x=x_j$.\n",
    "- When the current event deactivates, $x$ becomes zero, regardless of any previous events that may still be active.\n",
    "\n",
    "Three examples of special cases are shown in the figure below:\n",
    "\n",
    "<img src=\"figures/pacing-examples.svg\" width=\"400px\" />\n",
    "\n",
    "Three examples of pacing:\n",
    "_(Top)_ Three overlapping events.\n",
    "_(Middle)_ Event 1 deactivates just when event 2 activates, the pacing level never goes down to zero but jumps from $x_1$ to $x_2$.\n",
    "_(Bottom)_ An indefinitely re-occuring event with $d_i = p_i$: the pacing value stays at $x_i$ without interruption."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Why event-based pacing\n",
    "\n",
    "The two most common pacing methods used in cardiac cell simulation are (1) a periodic brief block current, or stimulus, and (2) a series of voltage steps.\n",
    "Both cases fit neatly into the event-based scheme.\n",
    "\n",
    "These protocols could also be implemented as part of the model - conceptually, they are outside the cell and therefore outside the model, but then again, so are the physical constants and external concentrations frequently included in cell models.\n",
    "\n",
    "However, a typical stimulus current is a 0.5ms pulse given once every 1000ms.\n",
    "For an efficient, variable step-size solver, this means there is a very large chance the solver will step over the stimulus completely.\n",
    "A simple way to remedy this is to limit the maximum step size to less than 0.5ms, but (1) this requires the solver to know about the protocol and (2) this degrades the efficiency of the variable step-size scheme, especially in the diastolic phases where step sizes can easily exceed 100ms.\n",
    "A much more efficient way to fix it is to tell the solver where the discontinuities are, so that it can ensure it doesn't step over them.\n",
    "This is what the event-based protocol achieves."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Time-series pacing\n",
    "\n",
    "In real experiments, an \"AP-clamp\" or \"data-clamp\" protocol is sometimes employed during which a pre-calculated time-series signal is applied to the cell, either as a current or a voltage (see e.g. [Bebarova, 2012](https://doi.org/10.4149/gpb_2012_016) ).\n",
    "In Myokit single-cell simulations, there are several ways to achieve something similar:\n",
    "\n",
    "- For piecewise constant signals, use the event-based pacing protocol. This is the most efficient method.\n",
    "- For sinusoidal or other time-dependent, continuous signals, add equations to the model.\n",
    "- For signals with a piecewise constant part and a more complicated part, split the simulation into two or more parts.\n",
    "- If there is no other option, use the _fixed-form pacing_ capabilities of the `Simulation` class.\n",
    " \n",
    "Fixed-form pacing allows you to specify a time series as a tuple $(\\{t_0, t_1, ... , t_{n}\\}, \\{x_0, x_1, ... x_{n}\\})$ where $t_i$ is a point in time and $x_i$ is the corresponding pacing value.\n",
    "The times array may contain duplicates, but must be strictly non-decreasing.\n",
    "At any time $t$, the value of the pacing variable $x$ is then defined by:\n",
    "\n",
    "- Searching for indices $i$ and $i+1$ such that $t_i \\leq t < t_{i+1}$.\n",
    "- If such an interval can be found, the value is calculated from values $x_i$ and $x_{i+1}$ using linear interpolation.\n",
    "- If the requested time is earlier than $t_0$, the value $x_0$ is used.\n",
    "- If the requested time is later than $t_n$, the value $x_n$ is used.\n",
    "\n",
    "Using bisection search, this can be implemented quickly.\n",
    "However, (1) the interpolation is still only an approximation of the \"true\" value of whatever process generated the time series.\n",
    "And (2), the interpolated function is not smooth: its derivatives may be discontinuities at any point $(t_i, x_i)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## No pacing\n",
    "\n",
    "In the absence of a pacing protocol (event based _or_ time series), any variables bound to a pacing label will have value $0.\n",
    "\n",
    "For simulations allowing multiple protocols (currently only the CVODES Simulation), this means that:\n",
    "- Variables bound to any pacing label defined when the simulation is created will be set to 0 if the protocol is subsequently removed.\n",
    "- If no protocol is set, i.e. the simulation is created with `protocol=None`, then the label `pace` is still used. This ensures backwards compatibility.\n",
    "- If the user explicitly sets an empty dict of pacing protocols \n",
    "\n",
    "For simulations using only the label `pace`:\n",
    "- If a variable is bound to `pace`, it will be set to 0 if no protocol is used (even if no protocol is given when the simulation is created).\n",
    "\n",
    "### Alternative\n",
    "\n",
    "An alternative idea was to have the pacing variable be undefined.\n",
    "In this case, any variables bound to the pacing variable would become unbound and retain their original, user-specified value.\n",
    "_However,_ this would allow variables to change from bound to unbound without (A) changes to the model or (B) a change in simulation engine.\n",
    "So by removing a protocol the result of running a simulation with `log=myokit.BOUND` would change."
   ]
  }
 ],
 "metadata": {
  "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.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}