Between Cardiac Electrical and Mechanical Activation Markers by Coupling Bidomain and Deformation Models

the material coordinates of the undeformed cardiac domain $\widehat{\varOmega }$ , by $\mathbf{x}=(x_1,x_2,x_3)^T$ the spatial coordinates of the deformed cardiac domain $\varOmega (t)$ at time t, by $\mathbf{F}(\mathbf{X},t)=\frac{\partial \mathbf{x}}{\partial {\mathbf {X}}}$ the deformation gradient. The cardiac tissue is modeled as a nonlinear hyperelastic material satisfying the steady-state force equilibrium equation

$\begin{aligned} \mathrm{Div}({{\mathbf {F}}} {\mathbf {S}})=\mathbf{0}, \qquad \mathbf{X} \in \widehat{\varOmega }. \end{aligned}$

(1)

The second Piola-Kirchoff stress tensor $\mathbf{S}=\mathbf{S}^{pas}+\mathbf{S}^{vol}+\mathbf{S}^{act}$ is the sum of passive, volumetric and active components. The passive and volumetric components are defined as

$S_{ij}^{pas,vol}=\frac{1}{2}\left( \frac{\partial W^{pas,vol}}{\partial E_{ij}}+\frac{\partial W^{pas,vol}}{\partial E_{ji}}\right) \quad i,j=1,2,3,$

where $\mathbf{E}=\frac{1}{2}(\mathbf{C}-\mathbf{I})$ is the Green-Lagrange strain tensor, $W^{pas}$ is an exponential strain energy function modeling the myocardium as an orthotropic hyperelastic material (see [14]), and $W^{vol}=K\left( J-1\right) ^2$ is a volume change penalization term, with K a positive bulk modulus and $J=det \mathbf{F}$ , added in order to model the myocardium as nearly incompressible.

(b) Mechanical Model of Active Tension. The active component $\mathbf{S}^{act}$ is given in terms of the active tension, that we assume developed along the myofiber direction only, ${\mathbf {S}}^{act}= T_a \frac{\widehat{\mathbf{a}}_l \otimes \widehat{\mathbf{a}}_l}{\widehat{\mathbf{a}}_l^T {\mathbf {C}}\, \widehat{\mathbf{a}}_l}$ , with $\widehat{\mathbf{a}}_l$ the fiber direction in the reference configuration and the biochemically generated active tension $T_a=T_a \left( Ca_i, \lambda , \frac{d\lambda }{dt} \right)$ is calcium, stretch ( $\lambda = \sqrt{\widehat{\mathbf{a}}_l^T{\mathbf {C}}\widehat{\mathbf{a}}_l}$ ), and stretch-rate ( $\frac{d\lambda }{dt}$ ) dependent, with the dynamics described by a system of ODEs proposed in Land et al. [21].

(c) Electrical Model of Cardiac Tissue: The Bidomain Model. We will use the following parabolic-elliptic formulation of the modified Bidomain model on the reference configuration $\widehat{\varOmega }\times (0,T)$ ,

$\begin{aligned} \left\{ \begin{array}{ll} c_m J~ \frac{\partial \widehat{v}}{\partial t} - {\mathrm{Div}} ( J~\mathbf{F}^{-1} D_i \mathbf{F}^{-T} \mathrm{Grad}\,(\widehat{v}+\widehat{u}_e)) + J~i_{ion}(\widehat{v},\widehat{\mathbf{w}},\widehat{\mathbf{c}},\lambda ) = 0 \\ - {\mathrm{Div}}(J~\mathbf{F}^{-1} D_i \mathbf{F}^{-T} \mathrm{Grad}\,\widehat{v}) - {\mathrm{Div}}(J~\mathbf{F}^{-1} (D_i+D_e) \mathbf{F}^{-T} \mathrm{Grad}\,\widehat{u}_e) = J~\widehat{i}_{app}^e \\ \frac{\partial \widehat{\mathbf{w}}}{\partial t} - {\mathbf {R}}_w(\widehat{v},\widehat{\mathbf{w}}) = 0, \displaystyle \frac{\partial \widehat{\mathbf{c}}}{\partial t} - {\mathbf {R}}_c(\widehat{v},\widehat{\mathbf{w}},\widehat{\mathbf{c}}) = 0.\\ \end{array} \right. \end{aligned}$

(2)

for the transmembrane potential $\widehat{v}$ extracellular potential $\widehat{u}_{e}$ gating and ionic concentrations variables $(\widehat{\mathbf{w}},\widehat{\mathbf{c}})$ . This system is completed by prescribing initial conditions, insulating boundary conditions, and applied current $~\widehat{i}_{app}^e$ , see [11] for further details. The conductivity tensors are given by

$\begin{aligned} \begin{array} {ll} D_{i,e}(\mathbf{x})&= \sigma ^{i,e}_l~ \mathbf{a}_l(\mathbf{x}) \mathbf{a}_l^T(\mathbf{x}) + \sigma ^{i,e}_t ~\mathbf{a}_t(\mathbf{x}) \mathbf{a}_t^T(\mathbf{x})+ \sigma ^{i,e}_n~ \mathbf{a}_n(\mathbf{x}) \mathbf{a}_n^T(\mathbf{x}), \end{array} \end{aligned}$

(3)

where $\sigma ^{i,e}_l,~\sigma ^{i,e}_t,~\sigma ^{i,e}_n~$ are the conductivity coefficients in the intra- and extracellular media measured along and across the fiber direction $\mathbf{a}_l, \mathbf{a}_t, \mathbf{a}_n$ .

(d) Ionic Membrane Model and Stretch-Activated Channel Current. The functions $I_{ion}(v,\mathbf{w},\mathbf{c},\lambda )$ ( $i_{ion}=\chi I_{ion}$ ), $R_w(v,\mathbf{w})$ and $R_c(v,\mathbf{w},\mathbf{c})$ in the Bidomain model (2) are given by the ionic membrane model by ten Tusscher et al. [28], available from the cellML depository (models.cellml.org/cellml). The ionic current is the sum $I_{ion}(v,\mathbf{w},\mathbf{c},\lambda ) = I_{ion}^m(v,\mathbf{w},\mathbf{c}) + I_{SAC}$ of the ionic term $I_{ion}^m(v,\mathbf{w},\mathbf{c})$ given by the ten Tusscher model and a stretch-activated channel current $I_{SAC}$ . This last current is modeled as in [22] as the sum of non-specific and specific currents $I_{SAC} = I_{SAC,n} + I_{Ko}$ .

3 Methods

Space Discretization. We discretize the cardiac domain with an hexahedral structured grid $T_{h_m}$ for the mechanical model (1) and $T_{h_e}$ for the electrical Bidomain model (2), where $T_{h_e}$ is a refinement of $T_{h_m}$ , i.e. $$h_m$$

is an integer multiple of $$h_e$$

. We then discretize all scalar and vector fields of both mechanical and electrical models by isoparametric $$Q_1$$

finite elements in space.

Time Discretization. The time discretization is performed by a semi-implicit splitting method, where the electrical and mechanical time steps could be different. At each time step,

(a)

given

at time

, solve the ODE system of the membrane model with a first order IMEX method to compute the new $w^{n+1}$ , $c^{n+1}$ .

(b)

given the calcium concentration $Ca_i^{n+1}$ , which is included in the concentration variables $c^{n+1}$ , solve the mechanical problems (1) and the active tension system to compute the new deformed coordinates $\mathbf{x}^{n+1}$ , providing the new deformation gradient tensor $\mathbf{F}_{n+1}$ .

(c)

given $w^{n+1}$ , $c^{n+1}$ , $\mathbf{F}_{n+1}$ and $J_{n+1}=\det (\mathbf{F}_{n+1})$ , solve the Bidomain system (2) with a first order IMEX method and compute the new electric potentials $v^{n+1},\ u_e^{n+1}$ with an operator splitting method.

We refer to [9] for more details.

Computational Kernels. Due to the employed space and time discretization strategies, at each time step, the main computational efforts consist of:

(i)

solving the nonlinear system deriving from the discretization of the mechanical problem by a Newton-GMRES-Algebraic Multigrid method, see [10];

(ii)

solving the two linear systems associated with the elliptic and parabolic Bidomain equations using the Conjugate Gradient method preconditioned by a Multilevel Additive Schwarz preconditioner studied in [23]. Our parallel simulations have been performed on a Linux cluster using the parallel library PETSc [6] from the Argonne National Laboratory.

Domain Geometry and Parameter Calibration. The cardiac domain $\widehat{\varOmega }=\varOmega (0)$ is the image of a cartesian slab using ellipsoidal coordinates, yielding a portion of truncated ellipsoid. The values of the orthotropic conductivity coefficients (see (3)) used in all the numerical tests are $\sigma ^i_l=3$ , $\sigma ^i_t=0.31525$ , $\sigma ^i_n=0.031525$ , $\sigma ^e_l=2$ , $\sigma ^e_t=1.3514$ , $\sigma ^i_n=0.6757$ , all expressed in $\mathrm{m}\varOmega ^{-1}\mathrm{cm}^{-1}$ . The parameters of the orthotropic strain energy function are given in [14]. The bulk modulus is $K=200\,\mathrm{kPa}$ .