technical-notes/1-4-opencl-simulation.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Technical note: OpenCL simulation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These technical notes describe the OpenCL simulation implemented in Myokit.\n",
    "This is designed for rectangular grids (or more complicated meshes) of _coupled cells,_ connected via a _diffusion current_.\n",
    "It is also possible to disable the diffusion current, and use it to perform single-cell simulations in parallel."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We assume each cell model contains a membrane potential variable, defined by an equation:\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{d V}{d t}=-\\frac{1}{C_{m}}(I_{\\mathit{ion}}+I_{\\mathit{stim}}+I_{\\mathit{diff}})\\label{eq:dv/dt}\n",
    "\\end{equation}\n",
    "\n",
    "where $I_\\mathit{diff}$ is a diffusion current\n",
    "\n",
    "\\begin{equation}\n",
    "I_{\\mathit{diff},i}=\\sum_{j\\in N_{i}}g_{ij}\\left(V_{i}-V_{j}\\right)\\label{eq:idiff}\n",
    "\\end{equation}\n",
    "\n",
    "and \n",
    "\n",
    "- $N_{i}$ is the set of neighbours of cell $i$;\n",
    "- $g_{ij}$ is the conductance between cell $i$ and its neighbour $j$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Specifying conductances\n",
    "\n",
    "As of version 1.33.0, Myokit provides three ways of specifying the conductances $g_{ij}$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rectangular, homogeneous conductances\n",
    "\n",
    "In this mode, we assume either a \"strand\" (in 1d) or a rectangular grid (2d) of cells, connected with a single conductance $g_x$ (1d) or a horizontal and a vertial conductance $(g_x, g_y)$ (2d)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Rectangular, heterogeneous conductances\n",
    "\n",
    "In this mode, we still assume a strand or a rectangular grid, but each conductance is set separately.\n",
    "\n",
    "In the 1d case, for $n$ cells we, this means we need a list of $n - 1$ conductances.\n",
    "\n",
    "In the 2d case, for $(m \\times n)$ cells, we need two matrices. One of $(m - 1 \\times n)$ conductances in the direction of the first coordinate, and one of $(m \\times n - 1)$ conductances in the direction of the second coordinate."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Arbitrary geometries\n",
    "\n",
    "In this mode, the simulation is initialised in 1d mode, so that each cell is identified by a single subscript $i$ or $j$.\n",
    "Next, the use specifies a list of connections"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heterogeneous parameters (scalar fields)\n",
    "\n",
    "The `SimulationOpenCL` allows _literal constants_ (constants that do not depend on other constants) to be set using a *scalar field*, i.e. to specify a value per cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Heterogeneous pacing\n",
    "\n",
    "In most simulations, only a few cells are paced.\n",
    "The `SimulationOpenCL` allows the paced cells to be specificied either explicitly (using a list of cell indices) or by specifying a rectangle of paced cells.\n",
    "\n",
    "All paced cells receive the same pacing signal.\n",
    "However, it is possible to set up more complex pacing e.g. by defining a stimulus current as $I_i = p_i * u(t)$ and setting $p_i$ to vary from cell to cell."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relationship to single-cell notation\n",
    "\n",
    "As before, we can write the system of ODEs for each cell as:\n",
    "\n",
    "\\begin{equation}\n",
    "\\dot{y}(t)=f\\left(y(t),u(t),t, p\\right)\n",
    "\\end{equation}\n",
    "\n",
    "but we now need to take into account that a cell i's external input $u_i$ depends on $y_i$, plus those of its neighbours:\n",
    "\n",
    "\\begin{equation}\n",
    "\\dot{y_i}(t)=f\\left(y_i(t),u_i(y, t), t, p_i\\right)\n",
    "\\end{equation}\n",
    "\n",
    "Here we write $f$ for each cell, instead of $f_i$, to indicate the function is the same (we use the same model for every cell).\n",
    "\n",
    "However, we write $u_i$ to indicate diffusion currents and pacing vary from cell to cell, and $p_i$ to indicate the possible use of scalar fields."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Relationship to monodomain model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The monodomain model (see [Leon & Horacek, 1991](https://doi.org/10.1016/0022-0736(91)90077-Y) or [Wikipedia](https://en.wikipedia.org/wiki/Monodomain_model )) can be used to simulate the electrophysiology of large 2 and 3d tissues.\n",
    "Instead of explicitly defining cells and connections between cells, this model assumes that each point in space can be described as having some cell properties, and some conductive properties.\n",
    "\n",
    "Under the monodomain model, the electrophysiology of each point in space is described by:\n",
    "\n",
    "\\begin{equation}\n",
    "\\chi \\left(\n",
    "c_m \\frac{\\partial V}{\\partial t}\n",
    "+ i_\\text{ion}\n",
    "+ i_\\text{stim}\n",
    "\\right)\n",
    "=\n",
    "\\frac{k}{k + 1}\n",
    "\\nabla \\cdot D_i \\nabla V\n",
    "\\end{equation}\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "where\n",
    "\n",
    "$V$ is the local membrane potential (e.g. in V),\n",
    "\n",
    "$i_\\text{ion}$ is the transmembrane current per unit area (e.g. in A/m^2),\n",
    "\n",
    "$i_\\text{stim}$ is the applied stimulus current per unit area (e.g. in A/m^2),\n",
    "\n",
    "$c_m$ is the membrane capacitance per unit area (e.g. in C/m^2, but usually given as $1\\mu F/\\textit{cm}^2$),\n",
    "\n",
    "$\\chi$ is the surface area of the membrane per unit of volume (e.g. in m^2/m^3),\n",
    "\n",
    "$D_i$ is a tensor representing intracellular conductivity (e.g. in S/m), and\n",
    "\n",
    "$k$ is the ratio between intra- and extracellular conductivity.\n",
    "\n",
    "Note that we have defined $i_\\text{stim}$ here to use the same sign convention as the other transmembrane currents, so that a negative current will excite the cell. This is different from many classical discussions, where the sign of this current is reversed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Rewriting, we find\n",
    "\n",
    "\\begin{align}\n",
    "\\frac{\\partial V}{\\partial t}\n",
    "&=\n",
    "-\\frac{1}{c_m}\n",
    "\\left(\n",
    "      i_\\text{ion}\n",
    "    + i_\\text{stim}\n",
    "    -\\frac{1}{\\chi}\n",
    "    \\frac{k}{k + 1}\n",
    "    \\nabla \\cdot D_i \\nabla V\n",
    "\\right) \\\\\n",
    "&= -\\frac{1}{c_m}(i_\\text{ion} + i_\\text{stim} + i_\\text{diff})\n",
    "\\end{align}\n",
    "\n",
    "so that\n",
    "\n",
    "\\begin{equation}\n",
    "i_\\text{diff}\n",
    "    = -\\frac{1}{\\chi} \\frac{k}{k + 1} \\nabla \\cdot D_i \\nabla V\n",
    "    = -\\nabla \\cdot \\tilde{D} \\nabla V\n",
    "\\end{equation}\n",
    "\n",
    "where \n",
    "\n",
    "\\begin{equation}\n",
    "    \\tilde{D} = \\frac{1}{\\chi}\\frac{k}{k + 1}D_i\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Assuming $\\tilde{D}$ does not vary spatially**, we can use [finite differences](https://en.wikipedia.org/wiki/Finite_difference#Higher-order_differences) to discretise this on a 1-dimensional regular grid, and obtain an expression for the diffusion current in node $i$\n",
    "\n",
    "\\begin{align}\n",
    "i_{\\text{diff}, i}\n",
    "    &= -\\nabla \\cdot \\tilde{D} \\nabla V \\\\\n",
    "    &= -\\tilde{D} \\frac{\\partial^{2}V}{\\partial x^{2}} \\\\\n",
    "    &= -\\tilde{D} \\left( \\frac{V_{i + 1} - 2 V_i + V_{i - 1}}{\\Delta x^{2}}\\right) \\\\\n",
    "    &=  \\tilde{D} \\frac{1}{\\Delta x^{2}}\\sum_{j \\in N}\\left(V_i-V_j\\right)\n",
    "\\end{align}\n",
    "\n",
    "where $\\Delta x$ is the spacing between subsequent points on the grid, and the sum is taken over the set of points neighbouring $i$.\n",
    "We can use this same equation for points at the edges, which is equivalent to assuming zero-flux boundary conditions.\n",
    "\n",
    "(Note that the $i$ and $j$ in the equations above denote a cell index, while the $i$ in $D_i$ is for _internal_.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this expression, we can write\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dV}{dt} = -\\frac{1}{c_m} (i_\\text{ion} + i_\\text{stim} + i_\\text{diff})\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For models with currents in $\\mu A / \\textit{cm}^2$ we can use this equation directly.\n",
    "To use unnormalised currents, we replace the specific capacitance with total capacitance, to find\n",
    "\n",
    "\\begin{equation}\n",
    "\\frac{dV}{dt} = -\\frac{1}{C_m} (I_\\text{ion} + I_\\text{stim} + I_\\text{diff})\n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Two-dimensions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In two dimensions, and assuming a diagonal conductivity tensor with values $\\tilde{D_x}$ and $\\tilde{D_y}$, we can follow a similar procedure on a grid with spacing $\\Delta x$ and $\\Delta y$ to find\n",
    "\n",
    "\\begin{align}\n",
    "i_\\text{diff}\n",
    "    = \\frac{\\tilde{D_x}}{\\Delta x^{2}} \\sum_{j \\in N_x}(V_i - V_j)\n",
    "    + \\frac{\\tilde{D_y}}{\\Delta y^{2}} \\sum_{j \\in N_y}(V_i - V_j) \\\\\n",
    "\\end{align}\n",
    "\n",
    "where $N_x$ and $N_y$ are the sets of neighbours in the directions corresponding to $\\tilde{D_x}$ and $\\tilde{D_y}$ respectively."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Converting parameters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Assuming a rectangular grid with spacing $\\Delta x$ in all directions and a spatially homogeneous conductance $\\tilde{D}$ or $D_{i}$, we can write\n",
    "\n",
    "\\begin{equation}\n",
    "i_\\text{diff}\n",
    "    = \\tilde{D} \\frac{1}{\\Delta x^{2}}\\sum_{j \\in N}\\left(V_{i}-V_{x}\\right)\n",
    "    = \\bar{g} \\sum_{j \\in N}\\left(V_{i}-V_{x}\\right)\n",
    "\\end{equation}\n",
    "\n",
    "to find\n",
    "\n",
    "\\begin{equation}\n",
    "\\bar{g} = \\frac{\\tilde{D}}{\\Delta x^2} = \\frac{1}{\\Delta x^2}\\frac{1}{\\chi}\\frac{k}{k + 1}D_i\n",
    "\\end{equation}\n",
    "\n",
    "where $\\bar{g}$ is cell-to-cell conductance per unit membrane area."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In unnormalised form, we obtain\n",
    "\n",
    "\\begin{equation}\n",
    "g = A_m \\frac{1}{\\Delta x^2}\\frac{1}{\\chi}\\frac{k}{k + 1}D_i\n",
    "\\end{equation}\n",
    "\n",
    "where $A_m$ is the surface area of the cell membrane encompassed by each \"node\" of our discretisation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Converting parameters: heterogeneity"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Please note, however, that $g_{ij}$ is defined *between* cells/nodes while $\\tilde{D}$ is defined *for* each node.\n",
    "This difference is unproblematic in spatially homogeneous simulations, but leads to some difficulties if hetereogeneity is introduced.\n",
    "\n",
    "For example, if the conductivity can vary spatially, then $\\nabla \\cdot \\tilde{D} \\nabla V \\neq \\tilde{D} \\frac{\\partial^{2}V}{\\partial x^{2}} $, and the equation for $g_{ij}$ becomes dependent on the spatial derivative of $\\tilde{D}$.\n",
    "\n",
    "Similarly, when using unnormalised models there is no way to choose an $A_m$ *per connection*, rather than per cell, unless we assume $A_m$ does not vary spatially."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solver\n",
    "\n",
    "The solver used in the `SimulationOpenCL` class is a simple forward-Euler method (although Rush-Larsen updates can be used for Hodgkin-Huxley gating equations, if desired).\n",
    "This has the advantage of being quite a _stable_ method, meaning it can take relatively large steps (larger than methods that take _slightly_ more complex explicit update steps, but not as large as _implicit_ methods).\n",
    "It also has a low computational cost _per step_, requiring only a single evaluation of the RHS.\n",
    "A major downside is the requirement to choose a good step size.\n",
    "\n",
    "There are some plans to build an adaptive method in Myokit, still based on forward Euler.\n",
    "However, it is worth noting that, as [Garcia-Molla et al.](https://doi.org/10.1016/j.compbiomed.2013.10.023) point out, multi-cell simulations are often of re-entrant arrhythmias in which there are always depolarising cells which require the smallest step size.\n",
    "In their work, they found fixed step-size methods often outperformed adaptive ones."
   ]
  }
 ],
 "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.9.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}