CEPS
Bidomain model



Equations

The famous bidomain model comes from the homogeneization of the quasi-static electric potential problem at the cell scale, where the intra- and extra-cellular potentials are defined. Assuming that the myocardium is locally constituted of a periodic pattern of cardiac cells, the bidomain problem reads:

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

Unknowns:

  • $v(\bm{x},t)$ is the transmembrane potential, in $\si{\milli\volt}$,
  • $u_\mathrm{e}(\bm{x},t)$ is the extracellular 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_\mathrm{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.
  • $\sigma_\mathrm{e}(\bm{x})$ the tensorial extracellular conductivity, in $\si{\milli\siemens\per\centi\metre}$. Similar to $\sigma_\mathrm{i}$

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.

The same numerical methods as for the monodomain problem are used to solve the bidomain problem.



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.

# Bidomain model:
# component index 0 represents the transmembrane potential
# component index 1 represents the extracellular potential
# ===============
# Output controls
# ===============
# Output name
output file name : ./results/bidomainAB2
# Output each 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 -0.1 -0.1 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.vtk
# 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 conductivityType
# * -1 for the region attribute means everywhere
# * g_l g_t g_n are the terms of the diagonalized conductivity tensor
# * conductivity types: 0 for intracellular, 1 for extracellular
conductivity : -1 2.4 0.8 0.8 0
conductivity : -1 2.4 0.8 0.8 1
# ===========
# Ionic model
# ===========
# Format: <ionicModelName> <regionAttribute> <component> <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
# Only the transmembrane potential is excited
stimulus: 200 15.0 0.1 1. -1 0
# Let's add a second stimulation, with a different spatial profile
# the last label is a reference to the function defined just below
stimulus : -1 15. 1.1 1 -1 0 stimFunc1
# definition of the spatial variation of stimulation
# format: <label> <C4|C5|CINF|CONSTANT> <center_x> <center_y> <center_z> <diameter> <amplitude>
# here we use a ~5mm wide bell-shaped C-infty profile, centerend around the
# given coordinates
stimulation profile : stimFunc1 CINF -1.48443 -19.2075 18.6753 1.
# ==================
# 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 : GMRES
@ GMRES
well, gmres



Outputs

Snapshots of the transmembrane voltage and extracellular potential 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.

These are the expected potential maps at t=7.5s when using the example file above.

References

FIXME bidomain Refs