CEPS  24.01
Cardiac ElectroPhysiology Simulator
Numerics

PDE Discretization

CEPS uses Finite Elements Methods (FEM) to discretize and solve Partial Differential Equations. We won't go much into details about FEMs in this page, so it is advised for newcomers to read introduction courses about FEM, in particular on how the variational formulation of a PDE is obtained.

Mesh

The computational domain is given as input in the form of a geometrical mesh, in which each region of the domain is split into basic shapes. CEPS current implementation only works with meshes formed with geometrical simplices, i.e.

  • segments for curves
  • triangles for surfaces
  • tetrahedra for volumes

CEPS will take into account attributes that are given on mesh nodes and mesh cells, which are identifiers associated to points when building the mesh in order to constitute specific regions (e.g. stimulation zones, tissue alteration, ionic models selection, etc.)

Possible mesh file formats are .mesh (generated by e.g. gmsh), .ele .node .edge sets of files generated by TetGen, .vtk legacy or .pvtu VTK files.

Finite Elements

Finite Elements are constituted of four ingredients:

  • Geometrical data with the nodes defining the shape of the element (a line, a square, a tetrahedron, etc.), and the jacobian $J$ of the transformation from a reference element to the element itself:

    \[ \begin{tikzpicture}[x=2cm,y=2cm] \usetikzlibrary{arrows} \node (p0) at (0,0){\textbullet}; \node (p1) at (1,0){\textbullet}; \node (p2) at (0,1){\textbullet}; \node (p3) at (0.5,0) {\textopenbullet}; \node (p4) at (0.5,0.5){\textopenbullet}; \node (p5) at (0,0.5) {\textopenbullet}; \node[below left] at (p0) {$(0,0)$}; \node[below right] at (p1) {$(1,0)$}; \node[above left] at (p2) {$(0,1)$}; \draw[-] (0,0)--(1,0)--(0,1)--cycle; \node (s0) at (2.0,0.2){\textbullet}; \node (s1) at (2.8,0.9){\textbullet}; \node (s2) at (3.1,0.1){\textbullet}; \node (s3) at (2.4,0.55){\textopenbullet}; \node (s4) at (2.95,0.5){\textopenbullet}; \node (s5) at (2.55,0.15){\textopenbullet}; \draw[-] (2,0.2)--(2.8,0.9)--(3.1,0.1)--cycle; \draw[-triangle 45,bend left] (1,0.6) to (1.8,0.6); \node at (1.4,0.85) {$J = \left(\partial_i F_j\right)_{i,j}$}; \end{tikzpicture} \tdplotsetmaincoords{75}{110} \quad\quad \begin{tikzpicture}[tdplot_main_coords,scale=2] \usetikzlibrary{arrows} \node (p0) at (0,0,0){\textbullet}; \node (p1) at (1,0,0){\textbullet}; \node (p2) at (0,1,0){\textbullet}; \node (p3) at (0,0,1){\textbullet}; \node (p4) at (0.5,0,0) {\textopenbullet}; \node (p6) at (0,0.5,0) {\textopenbullet}; \node (p7) at (0,0,0.5) {\textopenbullet}; \node (p5) at (0.33,0.33,0.33){\textopenbullet}; \node (p8) at (0,0.5,0.5) {\textopenbullet}; \node (p9) at (0.5,0,0.5) {\textopenbullet}; \node (p10) at (0.5,0.5,0) {\textopenbullet}; \node[left] at (p0) {$(0,0,0)$}; \node[below] at (p1) {$(1,0,0)$}; \node[below] at (p2) {$(0,1,0)$}; \node[above] at (p3) {$(0,0,1)$}; \draw[-] (0,0,0)--(1,0,0)--(0,1,0)--cycle; \draw[-] (0,0,0)--(0,0,1)--(0,1,0); \draw[-] (0,0,1)--(1,0,0); \draw[-triangle 45,bend left] (0,1,0.6) to (0,1.8,0.6); \node at (0,1.4,0.85) {$J = \left(\partial_i F_j\right)_{i,j}$}; \begin{scope}[xshift=2.5cm] \node (s0) at (0,0,0){\textbullet}; \node (s1) at (0.2,-0.3,0.4){\textbullet}; \node (s2) at (0,0.5,0.4){\textbullet}; \node (s3) at (-0.1,0.3,0.9){\textbullet}; \draw[-] (0,0,0) -- (0.2,-0.3,0.4) node [pos=0.5] {\textopenbullet} -- (0,0.5,0.4) node [pos=0.5] {\textopenbullet} -- (0,0,0) node [pos=0.5] {\textopenbullet}; \draw[-] (0,0,0) -- (-0.1,0.3,0.9) node [pos=0.5]{\textopenbullet} -- (0,0.5,0.4) node [pos=0.5] {\textopenbullet}; \draw[-] (-0.1,0.3,0.9) -- (0.2,-0.3,0.4) node [pos=0.5]{\textopenbullet}; \end{scope} \end{tikzpicture} \]

  • Nodal points on which are positionned the degrees of freedom where the unknowns are sought. The unknowns are usually the value of a function at this point, but it could also be its derivative.
  • Basis functions on which the variational formulation of the PDE are projected. These projections lead to the discrete formulation of the problem.
  • Quadrature formulae used to evaluate the integrals of the variational formulation. For now, the same quadrature rules are used for all elements.

Implemented finite elements

  • FEP1 (Lagrange): geometrical simplices (segments, triangles, tetrahedron) with continuous first order polynomial basis functions.
  • FEP2 (Lagrange): geometrical simplices (segments, triangles, tetrahedron) with continuous second order polynomial basis functions.

Boundary conditions

Boundary conditions can be given on different sets of elements for Neumann and Robin BC or nodes for Dirichlet BCs. Neumann and Robin BC are taken into account when computing the integrals on boundary elements by replacing the value of $\partial_n u$ by the appropriate value. Dirichlet BCs are imposed by zeroing the rows of the discretization matrix corresponding to Dirichlet nodes (with 1 on the diagonal). It is also possible to zero the column corresponding to these nodes to preserve the matrix symmetry (if any), but this comes to the expense of extra matrix copy operations that may slow down the execution more than necessary.



Numerical schemes for PDEs

We present here the numerical schemes that are used to solve the elements of the problem that involve spatial operators. For the discretization of the ODE that models the ionic exchanges, refer to the Numerics page.

The numerical schemes are expressed after finite element spatial discretization has been performed: $(M_{ij})_{ij}$ is the mass matrix of products of basis functions $\int \varphi_i\varphi_j$, and $(K_{ij}[\sigma])_{i,j}$ designates the stiffness terms corresponding to the discretization of the $\int\nabla\cdot(\sigma\nabla \varphi_i)\varphi_j$ terms.

Backward differentiation solvers (SBDF)

This is a multi-step method (for orders greater than 1). For a method of order $p$, we solve at each iteration:

\[ Mv^{n+1} + \beta_p \dfrac{\delta t}{A_\mathrm{m}C_\mathrm{m}}K[\sigma_i]v^{n+1} = \displaystyle\sum_{k=0}^{p-1} \alpha_{p,k}M\left[v^{n-k}+\dfrac{\delta t}{C_\mathrm{m}}\left(I_\mathrm{ion}^{n-k}+I_\mathrm{app}^{n-k}\right)\right], \]

with the coefficients

\[ \alpha = \begin{pmatrix} 1 &&& \\ \nicefrac{4}{3} & \nicefrac{-1}{3} &&\\ \nicefrac{18}{11} & \nicefrac{-9}{11} & \nicefrac{2}{11} & \\ \nicefrac{48}{25} & \nicefrac{-36}{25} & \nicefrac{16}{25} & \nicefrac{-3}{25} \end{pmatrix}\quad\quad \beta_p =\begin{pmatrix} 1 \\\nicefrac{2}{3}\\\nicefrac{6}{11}\\\nicefrac{12}{25} \end{pmatrix} \]

The first iterations where data at previous steps is not necessarily available are performed with lesser order $p$.

Runge-Kutta

Runge-Kutta solvers are currently not available.*

Several explicit Runge-Kutta schemes can be used (see numerics). The complicated part is to determine the evolution function $f^n_k$ from one sub step $v^n_k$ to another. We compute it as the solution of

\[ Mf^n_k = -K[\sigma_i]v^n_k + M\left[v^n + \delta t \displaystyle\sum_{j=0}{k-1} b_{jk} f^n_j\right], \]

where $b_{jk}$ are the matricial RK coefficients.

Crank-Nicolson

This method involves previous states at two steps:

\[ \begin{array}{l} Mv^{n+1} + \dfrac{1}{2}\;\dfrac{\delta t}{A_\mathrm{m}C_\mathrm{m}}K[\sigma_i]v^{n+1} = -\dfrac{\delta t}{2}K[\sigma_i]v^n + M\left[v^{n}+\dfrac{3}{2}\dfrac{\delta t}{C_\mathrm{m}}\left(I_\mathrm{ion}^{n}+I_\mathrm{app}^{n}\right) -\dfrac{1}{2}\dfrac{\delta t}{C_\mathrm{m}}\left(I_\mathrm{ion}^{n-1}+I_\mathrm{app}^{n-1}\right)\right],\\ \end{array} \]

A backward Euler step is performed at first iteration in order to compute $v^1$ (SBDF1).



ODE Solvers

Several classes of ODE solvers are implemented for ionic models but also for more generic functions. We consider in this section that we want to solve the ODE system

\[ y'=f(t,y) \]

where $y$ can be a vectorial function (e.g. ionic variables on all PDE degrees of freedom). Some solvers are written based on the assumption that $f$ can be split into a linear part $a$ and a non-linear part $b$:

\[ y'=a(t,y)y + b(t,y). \]

In the following $y^n$ represents the approximation of $y(t^n)$ and $\delta t$ the time step between $t^n$ and $t^{n+1}$.

Note: ODE time steps are constant troughout the simulation. No adaptive time stepping methods are currently implemented.

Runge-Kutta solvers

A bunch of explicit Runge-Kutta (RK) methods is coded in CEPS, from order 1 (which is Explicit Euler) to order 5. The evaluation of $y^{n+1}$ given $y^n$ is based on the computation of $n$ increments $k_i$ at various times between $t^n$ and $t^{n+1}$. Then, a linear combination of those increments is used to compute $y^{n+1}$. The algorithm is the following:

\[\begin{array}{l} \text{for }i\in\ \llbracket 1,n\rrbracket,\quad k_i = f\left(t+a_i\delta t, y^n+\delta t\sum_{j=1}^{n} b_{ij}k_j\right),\\ y^{n+1} = y^n + \delta t\sum_{i=1}^{n} c_ik_i, \end{array}\]

where the coefficients $a_i,\ b_{ij}$ and $c_i$ can be written using the Butcher formalism:

\[\begin{array}{c|c} A&B \\\hline &C \end{array}\]

The method is explicit if $B$ is a lower triangular matrix with 0-diagonal.

The implementation follows this formalism, so it is rather easy to add RK methods. Here are the keywords for the currently implemented methods:

Input file Keyword Order Method
RK 1 1 Explicit Euler
RK Heun2 2 Second order Heun method
RK 2 2 Usual Predictor-Corrector two steps method
RK 3 3 Kutta's third order method
RK Heun3 3 Heun's third order method
RK Ralston3 3 Ralston's third order method [1]
RK SSP3 3 Strong Stability Preserving Third order method [2]
RK 4 4 Classic RK 4th order
RK 38 4 Also a classic, the 3/8 method
RK Ralston4 4 Ralston's 4th order method, minimizes truncation error [1]
RK Nystrom5 5 5th order [3]
RK Fehlberg5 5 Fehlberg's 5th order from the RK-45 Fehlberg scheme. [4]


Backward differentiation solvers (Keyword: SBDF <order>)

These simple multistep solvers are based on Taylor expansions on several steps. For example, for a method of order 3, given values and derivatives at times $t^n,\ t^{n-1}$ and $t^{n-2}$, we have

\[\begin{array}{cccccccccccc}\displaystyle y^{n+1} &=& y^n &+& \delta t f^n &+& \frac{1}{2}\delta t^2(f^n)' &+& \frac{1}{6}\delta t^3(f^n)'' &+\mathcal{O}(\delta t^4),\\[+15pt] y^{n+1} &=& y^{n-1} &+& 2\delta tf^{n-1} &+& 2\delta t^2(f^{n-1})' &+&\frac{4}{3}\delta t^3(f^{n-1})''&+\mathcal{O}(\delta t^4),\\[+15pt] y^{n+1} &=& y^{n-2} &+& 2\delta tf^{n-2} &+& 2\delta t^2(f^{n-2})' &+&\frac{4}{3}\delta t^3(f^{n-2})''&+\mathcal{O}(\delta t^4). \end{array}\]

To achieve order 3, we find a linear combination of the three expansions that nullifies the terms in $\delta t^2$ an $\delta t^3$. This leads to

\[ y^{n+1} = \displaystyle \frac{18}{11} \left(y^n-\frac{1}{2}y^{n-1}+\frac{1}{9}y^{n-2} + \delta t\left(f^n-f^{n-1}+\frac{1}{3}f^{n-2}\right)\right). \]

SBDF solvers are implemented for order 1 to 4.


Semi-implicit Euler scheme (Keyword: FBE)

This scheme simply solves

\[ y^{n+1} = y^n + \delta t(a^ny^n{n+1} + b^n). \]

Note that there will trouble when $ a=\delta t^{-1}$.


Rush-Larsen solvers (Keyword: RL <order>)

Rush-Larsen (RL) solvers are exponential integrators using the decomposition of the evolution function in a linear and a non-linear part, and adapted to the ionic models. Please refer to [5] for explanations on the derivation of the scheme and a complete analysis.

The general formula is

\[ y^{n+1} = y^n + \delta t\,\varphi_1(\alpha_n\delta t)(\alpha_n y^n+\beta_n), \]

where

\[ \varphi_1(z) = \displaystyle\frac{\mathrm{e}^z-1}{z},\quad \varphi_1(0)=1, \]

and the coefficients $\alpha_n$ and $\beta_n$ are computed using data from multiple previous steps. If we note $a_n=a(t^n,y^n)$ and $b_n=b(t^n,y^n)$, coefficients are given depending on the order. RL schemes are implemented for order 1 to 4:

Order 1:

\[\begin{array}{rcl} \alpha_n &=& a_n,\\[+15pt] \beta_n &=& b_n. \end{array}\]

Order 2:

\[\begin{array}{rcl} \alpha_n &=& \frac{3}{2}a_n - \frac{1}{2}a_{n-1},\\[+15pt] \beta_n &=& \frac{3}{2}a_n - \frac{1}{2}a_{n-1}. \end{array}\]

Order 3:

\[\begin{array}{rcl} \alpha_n &=& \frac{23}{12}a_n - \frac{16}{12}a_{n-1} + \frac{5}{12}a_{n-2},\\[+15pt] \beta_n &=& \frac{23}{12}a_n - \frac{16}{12}a_{n-1} + \frac{5}{12}a_{n-2}\\[+15pt] &+& \frac{\delta t}{12} \left(a_nb_{n-1}-b_na_{n-1}\right). \end{array}\]

Order 4:

\[\begin{array}{rcl} \alpha_n &=& \frac{55}{24}a_n - \frac{59}{24}a_{n-1} + \frac{37}{24}a_{n-2} - \frac{9}{24}a_{n-3},\\[+15pt] \beta_n &=& \frac{55}{24}b_n - \frac{59}{24}b_{n-1} + \frac{37}{24}b_{n-2} - \frac{9}{24}b_{n-3}\\[+15pt] &+& \frac{\delta t}{12} \left(a_n(3b_{n-1}-b_{n-2}) - b_n(3a_{n-1}-a{n-2})\right). \end{array}\]


Exponential Adams-Bashforth ODE Solvers (Keyword: EAB <order>)

Another kind of exponential integrators... FIXME write here.

[6]

EAB schemes are implemented for order 1 to 4.


References

[1] Ralston, Anthony. "Runge-Kutta Methods with Minimum Error Bounds". Math. Compu.

[2] Dale E. Durran, “Numerical Methods for Fluid Dynamics”, Springer. Second Edition.

[3] Dormand, J.R., Prince, P.J. New Runge-Kutta algorithms for numerical simulation in dynamical astronomy. Celestial Mechanics

[4] Ernst Hairer, Syvert Nørsett, and Gerhard Wanner. Solving Ordinary Differential Equations I: Nonstiff Problems, second edition, Springer-Verlag, Berlin. ISBN 3-540-56670-8.

[5] Yves Coudière, Charlie Douanla Lontsi, and Charles Pierre. Rush-Larsen time-stepping methods of high order for stiff problems in cardiac electrophysiology. arXiv

[6] Yves Coudière, Charlie Douanla Lontsi, and Charles Pierre. Exponential Adams–Bashforth integrators for stiff ODEs, application to cardiac electrophysiology. Math Comput Simulat