Mathematical background

Mathematical background#

This page gives a self-contained introduction to the mathematical models that fenicsx-beat solves, and fixes the notation used consistently throughout the demos. It is a condensed version of the material in Sundnes, Lines, Cai, Nielsen, Mardal and Tveito, Computing the Electrical Activity in the Heart [SLC+07], which is the standard reference for this field and the source we follow below. If you want the full derivations, error analysis and numerical background, that book is the natural next stop.

1. From single cells to continuum tissue#

The heart muscle is made up of billions of individual, electrically excitable cells. Each cell is surrounded by a membrane that separates the intracellular space from the extracellular space, and actively maintains a potential difference across it. When a cell is stimulated, ion channels in the membrane open and close in a coordinated way, producing a rapid depolarization followed by a slower repolarization — together called an action potential. Neighbouring cells are electrically coupled through gap junctions, so an action potential triggered in one place propagates as a wave through the tissue.

Modelling every cell individually is computationally impossible for anything larger than a small tissue sample. Instead, we use a volume-averaging argument: at every point in space we define potentials and currents as averages over a small volume that contains many cells, but is still small compared to the size of the tissue. This is exactly the same idea used in continuum mechanics to go from molecules to continuous solids and fluids. Both the intracellular and extracellular spaces are then treated as continuous domains that coexist at every point in the tissue — this is the starting point of the bidomain model.

2. Notation used throughout this documentation#

The following symbols are used consistently in this background page, in the demos, and (as far as possible) map directly onto arguments in the beat API:

Symbol

Meaning

Corresponds to

\(u_i\), \(u_e\)

Intracellular / extracellular potential

—

\(v = u_i - u_e\)

Transmembrane potential

pde.state / solver.pde.state

\(s\)

Vector of gating / concentration state variables

ODE states array

\(M_i\), \(M_e\)

Intracellular / extracellular conductivity tensors

—

\(M\)

Effective (monodomain) conductivity tensor

M argument to beat.MonodomainModel

\(\chi\)

Membrane surface-to-volume ratio

chi argument to beat.stimulation.define_stimulus / beat.conductivities.define_conductivity_tensor

\(C_m\)

Membrane capacitance (per unit area)

C_m argument to beat.MonodomainModel

\(I_{ion}(v, s)\)

Total ionic current density through the membrane

Right-hand side of the cell model (fun)

\(I_{stim}\)

Externally applied stimulus current

beat.Stimulus / I_s

\(\theta\)

Splitting / time-stepping parameter, \(\theta \in [0, 1]\)

theta argument to beat.MonodomainSplittingSolver and parameters["theta"]

\(\theta = 1\) gives (first-order) Godunov splitting and a backward Euler PDE step, while \(\theta = 0.5\) gives (second-order) Strang splitting and a Crank–Nicolson PDE step; see Section 5 below.

3. The bidomain model#

At each point of the heart we now have two potentials, \(u_i\) and \(u_e\), and correspondingly two currents given by Ohm’s law,

\[ J_i = -M_i \nabla u_i, \qquad J_e = -M_e \nabla u_e, \]

where \(M_i\) and \(M_e\) are the (in general anisotropic) intracellular and extracellular conductivity tensors. Charge is conserved in each domain up to the current \(I_{ion}\) that crosses the membrane (measured per unit membrane area, so it is scaled by the surface-to-volume ratio \(\chi\) to become a current per unit tissue volume), giving

\[ -\nabla \cdot J_i = \chi C_m \frac{\partial v}{\partial t} + \chi I_{ion}(v, s), \qquad -\nabla \cdot J_e = -\chi C_m \frac{\partial v}{\partial t} - \chi I_{ion}(v, s), \]

where \(v = u_i - u_e\) is the transmembrane potential and \(C_m\) is the membrane capacitance. Eliminating \(u_i = u_e + v\) from these two equations gives the standard formulation of the bidomain model,

\[ \nabla \cdot (M_i \nabla v) + \nabla \cdot (M_i \nabla u_e) = \chi C_m \frac{\partial v}{\partial t} + \chi I_{ion}(v, s) - \chi I_{stim}, \]
\[ \nabla \cdot (M_i \nabla v) + \nabla \cdot \big((M_i + M_e) \nabla u_e\big) = 0, \]

together with an ODE system for the gating/concentration variables,

\[ \frac{\partial s}{\partial t} = f(v, s), \]

and insulating boundary conditions \(n \cdot (M_i \nabla v + M_i \nabla u_e) = 0\) and \(n \cdot (M_e \nabla u_e) = 0\) on the heart surface (if the heart is instead coupled to a surrounding conductive torso, these are replaced by continuity conditions on the potential and the current; see Chapter 2 of [SLC+07]). Note that we have added the stimulus current \(I_{stim}\) explicitly on the right-hand side — this is how an external pacing or defibrillation current enters the model, and it is treated identically to \(I_{ion}\) but with the opposite sign.

4. The monodomain model#

The bidomain model is a coupled system for two unknown fields, \(v\) and \(u_e\), and is comparatively expensive to solve. If we assume that the intracellular and extracellular conductivity tensors are proportional, \(M_e = \lambda M_i\) for some scalar \(\lambda\) (the equal anisotropy ratio assumption), then \(u_e\) can be eliminated algebraically and the bidomain model collapses to a single scalar reaction-diffusion equation for \(v\) alone, the monodomain model:

\[ \chi C_m \frac{\partial v}{\partial t} = \nabla \cdot (M \nabla v) - \chi I_{ion}(v, s) + \chi I_{stim}, \qquad M = \frac{\lambda}{1 + \lambda} M_i, \]

with the ODE system \(\partial s/\partial t = f(v, s)\) and the boundary condition \(n \cdot (M \nabla v) = 0\). The equal anisotropy assumption does not hold exactly in real tissue, so the monodomain model is an approximation — but it is far cheaper to solve and reproduces propagation patterns well enough that it is the default choice for most simulation studies, including all of the demos in this repository. beat.MonodomainModel implements exactly this equation. Note that the helper functions in beat.conductivities and beat.stimulation already divide the physical conductivities and stimulus amplitude by \(\chi\) before they reach MonodomainModel, so what the assembled weak form actually contains is

\[ C_m \frac{\partial v}{\partial t} = \nabla \cdot (M \nabla v) - I_{ion}(v, s) + I_{stim}, \]

with \(M\) and \(I_{stim}\) already expressed “per unit membrane area” rather than “per unit tissue volume”.

5. Ionic current models#

\(I_{ion}(v, s)\) describes the current flowing through the ion channels, pumps and exchangers embedded in the cell membrane, and it is here that most of the biophysical detail of a specific cell type enters the model. Three broad families are used in this repository’s demos:

  • Phenomenological models, such as the two-variable FitzHugh–Nagumo model used in the “Getting started” demo, are not derived from channel physiology but are simple dynamical systems tuned to reproduce the qualitative shape of an action potential (a fast upstroke, a plateau/recovery phase, and a refractory period). They are cheap to evaluate and useful for illustrating numerical methods and propagation patterns.

  • Detailed biophysical models in the tradition of Hodgkin and Huxley’s 1952 model of the squid giant axon describe the current through each ion channel as the product of a maximum conductance, a driving force \((v - v_{eq})\) relative to the channel’s Nernst equilibrium potential, and one or more dimensionless gating variables \(g \in [0, 1]\) (the fraction of open channels), each obeying a first-order ODE \(dg/dt = \alpha(v)(1 - g) - \beta(v) g\). The ten Tusscher–Panfilov and ToR-ORd models used in several demos (e.g. the Niederer benchmark, the pacing train, PVC, and the left/bi-ventricular ellipsoid demos) are modern, highly detailed descendants of this formalism, with tens of state variables describing individual sodium, potassium and calcium currents, calcium handling, and more.

  • Cell-model files for these detailed models are stored as CellML/.ode files under odes/ and translated into fast Python right-hand-side functions using gotranx, typically choosing a (generalized) Rush–Larsen scheme that treats the gating equations semi-analytically for numerical stability.

Regardless of which cell model is used, \(I_{ion}\) and \(s\) always play the same role in the coupled system above: \(I_{ion}(v, s)\) feeds back into the PDE for \(v\), while \(v\) itself feeds into the ODE system \(\partial s/\partial t = f(v, s)\) that advances the gating and concentration variables.

6. Solving the coupled system: operator splitting#

The full monodomain problem couples a (typically stiff) nonlinear ODE system at every point in space with a linear diffusion PDE. Instead of solving this monolithic system directly, fenicsx-beat uses operator splitting: at each time step, the ODE part and the PDE part are solved separately, using whichever numerical method suits each sub-problem best (an explicit or Rush–Larsen scheme for the stiff cell-model ODEs, and an implicit finite element solve for the diffusion PDE).

Writing the monodomain equation as \(\partial v/\partial t = (L_1 + L_2) v\) with

\[ L_1 v = -I_{ion}(v, s), \qquad L_2 v = \frac{1}{C_m}\nabla \cdot (M \nabla v) + \frac{1}{C_m} I_{stim}, \]

one time step \([t_n, t_n + \Delta t]\) of the general \(\theta\)-splitting scheme consists of three stages:

  1. ODE step. Solve \(\partial v/\partial t = -I_{ion}(v, s)\) together with \(\partial s/\partial t = f(v, s)\) from \(t_n\) to \(t_n + \theta \Delta t\).

  2. PDE step. Solve the linear diffusion equation \(C_m \partial v/\partial t = \nabla \cdot (M \nabla v) + I_{stim}\) from \(t_n\) to \(t_n + \Delta t\), using a \(\theta\)-rule (a weighted average of the diffusion operator at the old and new time levels) for the time discretization and a standard finite element method in space.

  3. Corrective ODE step. If \(\theta \ne 1\), solve the same ODE system as in step 1 again, from \(t_n + \theta \Delta t\) to \(t_n + \Delta t\), to bring \(v\) and \(s\) to the same point in time.

Choosing \(\theta = 1\) collapses step 3 and gives Godunov splitting (first-order accurate in time, and consistent with a first-order backward Euler PDE step); choosing \(\theta = 1/2\) gives Strang splitting combined with a Crank–Nicolson PDE step, which is second-order accurate overall. This is exactly what beat.MonodomainSplittingSolver implements: its theta argument selects between the two, and the pde and ode objects it is constructed from correspond to steps 2 and 1/3 above, respectively. The convergence and verification demos check this second-order behaviour numerically using the method of manufactured solutions.

7. Alternative backends for the ODE/PDE steps#

beat.MonodomainSplittingSolver only depends on its pde and ode objects implementing a small interface (step, to_dolfin/from_dolfin, …), so either half of the split can be swapped out independently. Two optional backends are available.

7.1 Fully implicit Runge–Kutta stepping (Irksome)#

The \(\theta\)-rule used in Section 6 gives at most second-order accuracy in time (at \(\theta = 1/2\)), and the explicit or Rush–Larsen schemes typically used for the ODE step in step 1/3 can struggle with very stiff cell models. beat.IrksomeMonodomainModel and beat.IrksomeODESolver are alternative implementations of the pde and ode roles in beat.MonodomainSplittingSolver, built on top of the Irksome library, that instead advance \(v\) (and, for the ODE step, \(s\)) with a fully implicit Runge–Kutta method specified by a Butcher tableau — for example irksome.BackwardEuler() (first order, matching Godunov splitting) or irksome.GaussLegendre(2) (fourth order). Because they implement the same interface as MonodomainModel and beat.odesolver.DolfinODESolver, an IrksomeMonodomainModel PDE step can be freely combined with a plain DolfinODESolver ODE step, or vice versa, inside the same splitting solver. IrksomeMonodomainModel can also be used entirely on its own, without any splitting at all, to solve a pure diffusion-plus-stimulus problem implicitly and with high-order accuracy in time — useful when there is no reaction term to split against, and hence no splitting error to bound the overall accuracy. beat.IrksomeMultiODESolver is the Irksome counterpart to beat.odesolver.DolfinMultiODESolver, for tissue with a different cell model in different regions (as in the multi-layer ventricle demos); it shares one Butcher tableau across all regions.

This comes at a cost: each implicit Runge–Kutta stage requires solving a (in general nonlinear) system over the whole mesh, so it is more expensive per time step than the explicit/\(\theta\)-rule schemes used elsewhere in this repository. irksome is therefore an optional dependency (the irksome extra, python3 -m pip install fenicsx-beat[irksome]), and the Irksome demo shows it in combination with a gotranx-generated cell model.

7.2 Cell-model evaluation via dolfinx-external-operator#

beat.ExternalOperatorODESolver and beat.ExternalOperatorMultiODESolver are drop-in alternatives to beat.odesolver.DolfinODESolver/DolfinMultiODESolver built on dolfinx-external-operator’s FEMExternalOperator, which lets an ODE right-hand side be evaluated as a plain NumPy or Numba array operation at every mesh point and spliced back into the dolfinx Function machinery, instead of DolfinODESolver’s own hand-rolled loop. Unlike the Irksome classes above, these keep the same fun(states, t, parameters, dt) cell-model convention already used throughout this repository unchanged — the same gotranx-generated or hand-written functions that work with DolfinODESolver work here too, with no separate UFL- or Irksome-specific variant needed. The numerical scheme (e.g. forward Euler, Rush–Larsen) and its order of accuracy are unchanged; what changes is only how it is evaluated and wired into the splitting solver, which mainly matters if you want to batch the cell-model evaluation on a backend dolfinx-external-operator supports (Numba, JAX, PyTorch) rather than beat’s own optional Numba path in single_cell.py/odesolver.py. dolfinx-external-operator is an optional dependency (the external_operator extra, python3 -m pip install fenicsx-beat[external_operator]), and the external-operator demo cross-checks it against DolfinODESolver on a real gotranx-generated cell model and reports the current performance trade-off honestly.

7.3 Removing the splitting error entirely: monolithic implicit coupling#

Sections 7.1 and 7.2 both still fit inside the operator-splitting scheme of Section 6 — they change how the PDE or ODE half-step is computed, not the fact that there are two half-steps with a splitting error between them. Because FEMExternalOperator can supply a derivative as well as a value, it can instead be used to solve the PDE and the cell-model ODEs together, in one fully implicit Newton system per time step, with \(v\) and every gating/concentration variable as unknowns simultaneously — removing the operator-splitting error altogether (only the usual temporal and spatial discretization error remains). This needs a cell model that can supply its own Jacobian; gotranx’s JAX code generation backend combined with jax.jacfwd/jax.vmap gives this for free, with no hand-derived derivatives, even through the Conditional-based branching common in cell models.

This is not currently packaged as a reusable beat class the way the classes above are — FEMExternalOperator does not (yet) support operands built from a mixed-space Function, so a monolithic coupling has to be built from one separate scalar Function per state plus an \(N \times N\) block Jacobian, which does not scale to large state counts as cleanly as the flat (num_states, num_points) arrays used everywhere else in this repository. The monolithic coupling demo is a validated prototype of the pattern (checked against an independent, non-FEM reference) together with a quantitative comparison against MonodomainSplittingSolver, showing a several-times reduction in error at practical time steps — worth knowing the pattern exists even though it is not (yet) a drop-in solver class.

8. Where to go from here#

  • The demos index lists every demo in this repository, grouped by topic.

  • The FitzHugh–Nagumo demo is the best starting point: it solves the monodomain model with the simplest possible ionic current model, both as a single-cell ODE and as a full 2D tissue simulation, and introduces the beat API step by step.

  • The diffusion demo shows the monodomain model with the ionic current turned off entirely, i.e. pure diffusion driven by a stimulus.

  • The convergence and verification demos connect the numerical scheme described in Section 6 back to measured convergence rates.

  • The remaining demos build up towards realistic simulations: conduction velocity and ECG estimation on a simple slab, pacing trains and premature ventricular complexes on a 1D cable, and full biophysically detailed simulations on left- and bi-ventricular geometries, including 12-lead ECG recovery.

  • The Irksome demo shows the fully implicit, high-order time stepping described in Section 7.1, the external-operator demo shows the drop-in ODE-step backend described in Section 7.2, and the monolithic coupling demo shows the fully implicit, no-splitting-error prototype described in Section 7.3.

  • For the full derivation of the bidomain and monodomain models, the physiology of the cell membrane, the Hodgkin–Huxley and Nernst–Planck formalisms, and a rigorous treatment of operator splitting and its order of accuracy, see Chapters 1–3 of [SLC+07].

[SLC+07] (1,2,3)

Joakim Sundnes, Glenn Terje Lines, Xing Cai, Bjørn Frederik Nielsen, Kent-Andre Mardal, and Aslak Tveito. Computing the electrical activity in the heart. Volume 1. Springer Science & Business Media, 2007.