Cardiac Electrophysiology Simulation using svFSI

The electrical activity in a heart is initiated at the sinoatrial node, which is located at the right atrium of the heart. A small cluster of pacemaker cells periodically emits electrical signals into the heart’s conduction system. These signals first travel at a high speed through a network of special cells that allows the signals to reach the whole heart rapidly. Finally, they travel into the ventricular muscles, and cause the entire heart to contract almost simultaneously. The electrical signals travel between cells through a process called depolarization, during which voltage-sensitive protein channels on the cell membrane open to allow positively charged ions to move in and out of the cell, leading to ionic currents.

Mathematical Model of Cardiac Electrophysiology

Propagation model

The propagation of electrical signal in the heart is governed by a reaction-diffusion like equation,

$$\frac{\partial V}{\partial t} + \frac{I_{ion} - I_{app}(t)}{C_m} = \nabla \cdot\left( \mathbf{D}\nabla V \right)$$ $$\mathrm{in~} \Omega^E\times(0,T] \nonumber$$ $$(\mathbf{D}\nabla V) \cdot \mathbf{N}=0 \mathrm{~on~} \partial\Omega^E\times(0,T]$$

$V$ is the local transmembrane potential. $C_m$ is the local membrane capacitance per unit area. $I_{ion}$ and $I_{app}$ are the ionic current flux (current per unit area) and applied external current flux, respectively. Here, $\mathbf{D}$ dictates the propagation velocity of the electrical signal and has the similar physical meaning as the diffusivity. It is calculated as

$$\mathbf{D} =\frac{\sigma}{\chi C_m}$$

$\sigma$ is the conductivity tensor and has the unit $S/m$. $\chi$ is the ratio of membrane area over tissue volume. Then, $\mathbf{D}$ is a tensor with the dimension $Length^2/Time$.

The above equation is the mono-domain description of the cardiac electrophysiology, i.e. we don’t solve the intra- and extra-cellular electrical signal propagation separately. The mono-domain and multi-domain conductivities are connected through the following relation,

$$\sigma = \frac{\sigma_i\sigma_e}{\sigma_i + \sigma_e}$$

where $\sigma_i$ and $\sigma_e$ are the intra- and extra-cellular conductivity tensors [1]. It is commonly assumed that the conductivity is transversely isotropic,

$$\mathbf{\sigma} = \sigma_f \mathbf{f}\otimes \mathbf{f} + \sigma_s (\mathbf{I}-\mathbf{f}\otimes \mathbf{f})$$

where $\sigma_f$ and $\sigma_s$ are the conductivities along the fiber direction and in the transverse plane. $\mathbf{f}$ is the fiber direction vector.

In svFSI, we directly specify $\mathbf{D}$ in the input file. We choose a slightly different form to enforce the transverse isotropy of the conductivity tensor,

$$\mathbf{D} = D_{iso}\mathbf{I} + \sum_{n=1}^{nsd}D_{ani,n}\mathbf{fN}_n\otimes\mathbf{fN}_n$$

Here, $nsd$ is the dimension, and $\mathbf{fN_n}$ is the local orthonormal coordinate system built by fiber direction and sheet direction. To connect with the previous expression, we have $D_f=D_{iso}+D_{ani}$ and $D_s=D_{iso}$.

Cellular activation models

Depending on how the depolarization and repolarization within a single cardiac myocyte is described, the electrophysiology models can be roughly divided into two categories: biophysics-based ionic models (such as the ten Tusscher-Panfilov (TTP) model[2][3]), and phenomenological models (such as the Aliev-Panfilov (AP), Fitzhugh-Nagumo (FN) models[4]).

Biophysics-based ionic models

Between cardiac myocytes, the propagation of the electrical signal is enabled by the transmembrane motion of different ions, such as $K^+$, $Na^+$ and $Ca^{2+}$ during depolarization and repolarization of the myocytes. Cellular biophysics-based activation models are designed to capture these ionic movements. One of the popular biophysics-based ionic models, the TTP model [3] is implemented in svFSI. Here, we will briefly illustrate the ionic currents involved in the TTP model.

In the TTP model, the ionic current is expressed as

$$I_{ion} = I_{Na} + I_{Kl} + I_{to} + I_{Kr} + I_{Ks} + I_{CaL} + I_{NaCa} + I_{NaK} + I_{pCa} + I_{pK} + I_{bCa} + I_{bNa}$$

The governing equations for each current can be found in [3].

Since, intracellular calcium concentration is the driving factor behind excitation-contraction coupling, we focus on the calcium dynamics here. In the TTP model,the following calcium concentrations are included:

• $Ca_{itotal}$: total (free+buffered) cytoplasmic $Ca^{2+}$ concentration.
• $Ca_{SRtotal}$: total SR $Ca^{2+}$ concentration.
• $Ca_{SStotal}$: total dyadic space $Ca^{2+}$ concentration.
• $Ca_i$: free cytoplasmic $Ca^{2+}$ concentration.
• $Ca_{SR}$: free SR $Ca^{2+}$ concentration.
• $Ca_{SS}$: free dyadic space $Ca^{2+}$ concentration.

Sarcoplasmic reticulum (SR) is a membrane structure within the muscle cell, whose main function is to store $Ca^{2+}$. Dyadic space is the region bounded by the T-tubule and SR. $Ca^{2+}$ ions are considered buffered if they are bound to negatively charged proteins (called buffers). Otherwise they are considered free. The action potential transmitted through the gap junction cause the current myocyte to depolarize. The depolarization opens the L-type (long-lasting) Ca channel located on the surface membrane. A small amount of $Ca^{2+}$ enters the myocyte due to potential, leading to a sharp increase of $Ca^{2+}$ in the dyadic space, which is a small region. This increase triggers the SR to release a large amount of $Ca^{2+}$ (calcium-induced calcium release) to enable the excitation-contraction. During diastole, calcium is removed from the cytoplasm through two ways. Ca is pumped (1) back into the SR and (2) out of the cell, mainly by the sodium-calcium exchange (NCX).

For the TTP model in svFSI, the following units have to be used: time in [ms], length [mm], amount of substance [mmol], voltage [mV]. Mass can be in [g].

Phenomenological models

Phenomenological models are derived based on some observations of the full ion model [4]. Instead of following the transmembrane ionic currents, they use an oscillation system with a fast ($V$) and a slow ($r$) variable to mimic the behaviors of the action potentials. The oscillators, without considering the diffusion term, are modeled as

$$\frac{\mathrm{d}V}{\mathrm{d}t}=f^{V}(V,r)$$ $$\frac{\mathrm{d}r}{\mathrm{d}t}=f^{r}(V,r)$$

Note that this set of equations describes the electrophysiology in a single cardiac myocyte, and the choice of $f^V$ and $f^r$ will determine if this is a pacemaker cell or a non-pacemaker cell.

The FitzHugh-Nagumo (FN) model describes the pacemaker cells with

$$f^V = c[V(V-\alpha)(1-V)-r]$$ $$f^r = V-br+a$$

The Aliev-Panfilov (AP) model describes the non-pacemaker cells with

$$f^V = cV(V-\alpha)(1-V)-Vr$$

$$f^r = \left( \gamma+\frac{\mu_1r}{\mu_2+V}\right)[-r-cV(V-b-1)]$$

Note that $V$ and $t$ are non-dimensional values here. The following equations are used to recover the physiological action potential and time:

$$\mathrm{FitzHugh-Nagumo model}: V^{fhn}=(65V-35)mV; ~~ t^{fhn} = (220t) ms$$

$$\mathrm{Aliev-Panfilov model}: V^{ap}=(100V-80)mV; ~~ t^{ap} = (12.9t) ms$$

Another class of phenomenological models exists that include additional variables to account for intra-cellular calcium kinetics. One such model is the Bueno-Orovio (BO) model [5] that can serve as a trade-off between the complex TTP model and the simplified phenomenological AP model.

Activation models available in svFSI

The following table provides a summary of all the available electrophysiology models in svFSI.

Available cardiac electrophysiology models.
Electrophysiology Model Input Keyword
Aliev-Panfilov model[4] “ap”, “aliev-panfilov”
Fitzhugh-Nagumo model[4] “fn”, “fitzhugh-nagumo”
Bueno-Orovio-Cherry-Fenton model[5] “bo”, “bueno-orovio”
tenTusscher-Panfilov model[3] “ttp”, “tentusscher-panfilov”

Example - Activation of a Myocardial Block

In this section, we will provide a example of cardiac electrophysiology modeling. This example is a widely used benchmark test, activation of a block of cardiac tissue [1]. More examples are available on GitHub.

The computational domain is a rectangular block, with dimensions of $3\times 7\times 20 mm$. A small cluster of cells located at one corner of the block (marked as S) will recive an initial stimulus, and the TTP model is used to model the activation of the rest of the domain. The entire simulation set-up along with results can be found on GitHub.

Here is the breakdown of the input file for this study.

# File svFSI.inp
#----------------------------------------------------------------
# General simulation parameters

Continue previous simulation: f
Number of spatial dimensions: 3
Number of time steps: 100000
Time step size: 0.01
Spectral radius of infinite time step: 0.50
Searched file name to trigger stop: STOP_SIM

Save results to VTK format: 1
Name prefix of saved VTK files: result
Increment in saving VTK files: 100
Start saving after time step: 1

Increment in saving restart files: 100
Convert BIN to VTK format: 0

Verbose: 1
Warning: 1
Debug: 0

#----------------------------------------------------------------
# Mesh data
Mesh file path):   mesh/h0.1/mesh-complete.mesh.vtu
Face file path: mesh/h0.1/mesh-surfaces/X0.vtp
}
Face file path: mesh/h0.1/mesh-surfaces/X1.vtp
}
Face file path: mesh/h0.1/mesh-surfaces/Y0.vtp
}
Face file path: mesh/h0.1/mesh-surfaces/Y1.vtp
}
Face file path: mesh/h0.1/mesh-surfaces/Z0.vtp
}
Face file path: mesh/h0.1/mesh-surfaces/Z1.vtp
}

Fiber direction: (1, 0, 0)

# Here we also need to supply a domain information to separate
# pacemaker cells and non-pacemaker cells. domain_info.dat
# is a plaintext file store the domain ID of each element.
# id == 1: non-pacemaker cells
# id == 2: pacemaker cells (recivies stimulus)
Domain file path: mesh/h0.1/domain_info.dat
}

#----------------------------------------------------------------
# Equations
Coupled: 1
Min iterations: 1
Max iterations: 5
Tolerance: 1e-6

# Here domain id == 1 are non-pacemaker cells
Domain: 1 {
Electrophysiology model: TTP
Conductivity (iso): 0.012571
Conductivity (ani): 0.082715
ODE solver: RK
}

# Here domain id == 2 are pacemaker cells
Domain: 2 {
Electrophysiology model: TTP
Conductivity (iso): 0.012571
Conductivity (ani): 0.082715
# Stimulus to initiate the depolarization process
Stimulus: Istim {
Amplitude: -35.714
Start time: 0.0
Duration: 2.0
Cycle length: 10000.0
}
ODE solver: RK
}

Output: Spatial {
Action_potential: t
}

LS type: GMRES
{
Max iterations:      100
Tolerance:           1D-6
Krylov space dimension: 50
}
}