CEPS
Monodomain model



Equations

The monodomain equation is a variant of the bidomain equation where the intra- and extracellular conductivities are assumed to be proportional:

\[ \sisetup{math-micro=\text{µ},text-micro=µ} \sigma_i = \lambda\sigma_e. \]

On the computation domain $\Omega$, the problem reads:

\[ \left\{\begin{array}{ll} A_\mathrm{m}(C_\mathrm{m} \partial_t v + I_\mathrm{ion}(v,w) + I_\mathrm{app}) - \dfrac{1}{1+\lambda}\nabla\cdot(\sigma_i\nabla v) = 0 & \text{on }\Omega,\\ \partial_t w + g(v,w) = 0 &\text{on }\Omega,\\ \sigma_i\partial_{\bm{n}} v = 0&\text{on }\partial\Omega,\\ + \text{ initial conditions.} \end{array}\right. \]

Unknowns:

  • $v(\bm{x},t)$ is the transmembrane potential, in $\si{\milli\volt}$,
  • $w(\bm{x},t)$ is a vector of ionic variables describing the state of the membrane.

Tissue properties:

  • $A_\mathrm{m}$ is the surface to volume ratio of cell membranes, in $\si{\per\centi\metre}$,
  • $C_\mathrm{m}$ is the cell membrane surfacic capacitance, in $\si{\micro\farad\per\centi\metre\squared}$,
  • $\sigma_i(\bm{x})$ the tensorial intracellular conductivity, in $\si{\milli\siemens\per\centi\metre}$. Usually, the tensor is the same everywhere when expressed in the local frame defined by fiber orientation. However, alterations of the tissue conductivity can be introduced by using volume fractions.

Fonctions:

  • $I_\mathrm{ion}$ is the current due to ionic gates, in $\si{\micro\ampere\per\centi\metre\squared}$,
  • $I_\mathrm{app}(\bm{x},t)$ is the user-defined applied current used to stimulate a specific region of the domain,
  • $g$ is the evolution function of ionic variables.



Numerical schemes for the PDE

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: $(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

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).



How to use

Here is an example of input file, which is available in the examples directory. See the input file page for more details on options.

#-----------------------------------------------------------
# Monodomain model:
# component index 0 represents the transmembrane potential
#-----------------
# Output controls
#-----------------
# Output name
output file name : ./results/monodomainAB2
# Output every 0.5ms
snapshot time : 0.5
# Uncomment to plot currents
#save ionic current : true
#save stimulation current : true
# Data used to generate the activation time map. Values are
# - AP detection threshold
# - Minimum voltage used for peak detection
# - AP duration percentage
activation time data : 0.2 0.2 0.75
# Probe points: values of the unknowns at those points*
# along time will be written in a separate file
# (*at closest point in given mesh to be exact)
probe points : 2.68655 -15.8401 20.9223 , 1.34332 -16.1237 19.8474
#-----------------
# Simulation time
#-----------------
# Times are given in milliseconds
PDE start time : 0.0
PDE end time : 25
#----------
# Geometry
#----------
# Meshes
2d mesh : ./meshes/atriaSurfacic
# The given mesh has units of cm, which is also CEPS length unit.
# So there is no need for scaling.
# geometry scale : 1
# Partitioning method can be
# - ParmetisElem
# - ParmetisElemGeom
# - PTScotchNode
partitioning method : PTScotchNode
#--------
# Tissue
#--------
# Surface to volume ratio, in cm-1
Am : 1.0
# Surfacic membrane capacitance, in muF per cm2
Cm : 1.0
# Conductivity, in mS per cm
# Format: regionAttribute g_l g_t g_n
# -1 for the region means everywhere
# g_l g_t g_n are the terms of the diagonalized conductivity tensor
conductivity : -1 1.2 0.4 0.4
#-------------
# Ionic model
#-------------
# Format: <ionicModelName> <regionAttribute> <componentIndex> <odeSolver>
# Here we use the Mitchell Schaeffer model everywhere (-1)
# for transmembrane potential (component 0) and solve the ODEs
# using the exponential Adams-Bashforth scheme of order 2.
ionic model : MS -1 0 EAB 2
#-------------
# Stimulation
#-------------
# Format: <attr> <amp> <tstart> <duration> <period> <component>
# Here we stimulate on points with attribute 200 (field "stimRegion" in vtk mesh file)
# Stimulations starts at 0.1ms and lasts 1ms. It happens once (period -1)
# Amplitude is in muA per cm2
stimulus: 200 15. 0.1 1 -1 0
#--------------------
# PDE solver options
#--------------------
finite elements : P1
PDE solver : SBDF 2
PDE time step : 0.05
# Linear solver options
absolute tolerance : 1.E-12
relative tolerance : 1.E-12
solver type : CG
@ CG
conjugate gradient



Outputs

Snapshots of the transmembrane voltage are written by the application. Ionic and stimulation currents can be written as well if needed. In addition, activation maps can be generated (see parameters above), as well as traces in time of the transmembrane voltage in a separate text file.

This is the expected activation time maps when using the example file above.



References

FIXME refs monodomain