An orthotropic electro-viscoelastic model for the heart with stress-assisted diffusion

TL;DR

This paper introduces a comprehensive orthotropic electro-viscoelastic model for cardiac tissue that integrates nonlinear elasticity, stress-dependent electrical conductivities, and electrophysiological dynamics, validated through numerical simulations.

Contribution

It presents a novel coupled electromechanical model with stress-assisted diffusion and a mixed-primal finite element scheme for cardiac tissue simulation.

Findings

Model accurately describes cardiac action potential

Stress influences electrical conductivities in the model

Numerical validation confirms model's robustness

Abstract

We propose and analyse the properties of a new class of models for the electromechanics of cardiac tissue. The set of governing equations consists of nonlinear elasticity using a viscoelastic and orthotropic exponential constitutive law (this is so for both active stress and active strain formulations of active mechanics) coupled with a four-variable phenomenological model for human cardiac cell electrophysiology, which produces an accurate description of the action potential. The conductivities in the model of electric propagation are modified according to stress, inducing an additional degree of nonlinearity and anisotropy in the coupling mechanisms; and the activation model assumes a simplified stretch-calcium interaction generating active tension or active strain. The influence of the new terms in the electromechanical model is evaluated through a sensitivity analysis, and we provide numerical validation through a set of computational tests using a novel mixed-primal finite element scheme.

Tables

Table 1. (a) Hyperelasticity variables

DoF	$h$	${\Vert 𝚷-{𝚷}_h\Vert }_{0,\Omega }$	rate	${\Vert 𝒖-{𝒖}_h\Vert }_{1,\Omega }$	rate	${\Vert p-p_h\Vert }_{0,\Omega }$	rate
$l=0$
77	0.7071	43.252	–	0.0576	–	30.161	–
253	0.3536	27.137	0.6725	0.0342	0.6345	19.030	0.6647
917	0.1768	12.535	1.1140	0.0216	0.7615	9.2110	1.0471
3493	0.0884	6.2636	1.0012	0.0118	0.8751	4.8012	0.9401
13637	0.0442	1.9169	1.1727	0.0071	0.9516	1.9631	1.3817
53893	0.0221	0.9841	0.9907	0.0042	0.9737	0.9206	0.9858
$l=1$
221	0.7071	19.481	–	0.0146	–	6.0355	–
789	0.3536	7.9032	1.3034	0.0037	1.7593	1.5809	1.4581
2981	0.1768	2.6409	1.8079	0.0011	1.7809	0.4120	1.7269
11589	0.0884	0.7277	1.9033	4.11E-4	1.8065	0.1353	1.8813
45701	0.0442	0.2063	1.9182	1.09E-4	1.9330	0.0382	1.9602
181509	0.0221	0.0569	1.9466	3.12E-5	1.9522	0.0094	1.9571

Table 2. Table 4.2: Model parameters for the electro-viscoelastic model ( 2.8 ), ( 2.11 ), ( 2.16 ), ( 2.14 ). Values are taken from cherubini17 ; gao ; rossi14 ; bueno08 , and the transmembrane potential v 𝑣 v is dimensionless.

Viscoelasticity constants

$a=$

0.236

[N/cm2]

$a_f=$

1.160

[N/cm2]

$a_s=$

3.724

[N/cm2]

$a_{fs}=$

4.010

[N/cm2]

$b=$

10.81

[–]

$b_f=$

14.15

[–]

$b_s=$

5.165

[–]

$b_{fs}=$

11.60

[–]

$p_0=$

0.1

[N/cm2]

$\beta =$

10

[ms]

$\delta =$

22.6

[N/cm2 ms]

${\zeta }_{\text{stab}}=$

0.25

[–]

${\eta }_a=$

0.001

[N/cm2]

${\eta }_b=$

0.01

[N/cm2]

${\kappa }_{sn}=$

0.6

[–]

${\kappa }_{nn}=$

0.03

[–]

${\rho }_0=$

0.001

[N/cm2]

Electrophysiology constants

$v_0=$

0

[–]

$v_v=$

1.55

[–]

$v_2^{-}=$

0.03

[–]

$v_{so}=$

0.65

[–]

$v_3=$

0.908

[–]

${\theta }_1=$

0.3

[–]

${\theta }_1^{-}=$

0.006

[–]

${\theta }_o=$

0.006

[–]

${\theta }_2=$

0.13

[–]

$k_2^{-}=$

65

[–]

$k_3=$

2.099

[–]

$k_{so}=$

2.045

[–]

$r_{2,\infty }^{\ast }=$

0.94

[–]

${\tau }_{2,\infty }=$

0.07

[–]

${\tau }_{1,1}^{-}=$

60

[–]

${\tau }_{1,2}^{-}=$

1150

[–]

${\tau }_{2,1}^{-}=$

60

[–]

${\tau }_{2,2}^{-}=$

15

[–]

${\tau }_{fi}=$

0.11

[–]

${\tau }_{o,1}=$

30.02

[–]

${\tau }_{o,2}=$

0.996

[–]

${\tau }_{so,1}=$

2.046

[–]

${\tau }_{so,2}=$

0.65

[–]

${\tau }_{3,1}=$

2.734

[–]

${\tau }_{3,2}=$

16

[–]

${\tau }_{si}=$

1.888

[–]

${\tau }_1^{+}=$

1.451

[–]

${\tau }_2^{+}=$

200

[–]

$\chi =$

1

[–]

Activation and excitation-contraction coupling constants

$D_0=$

1.171

[mm2/s]

$D_1=$

0.9

[mm2/s]

$D_2=$

0.01

[mm2/s]

$K_0=$

5

[–]

$K_1=$

-0.015

[–]

$K_2=$

-0.15

[–]

${\alpha }_1=$

10

[–]

${\alpha }_2=$

0.5

[–]

Table 3. Table 4.3: Convergence of conduction velocity with respect to temporal and spatial discretisation.

Convergence of Conduction Velocity, mm/ms
DoF	h(mm)	$\Delta t=$	$0.3\text{ms}$	$0.1\text{ms}$	$0.05\text{ms}$	$0.01\text{ms}$
27038	0.3817		0.1130	0.1032	0.1015	0.0994
108576	0.1909		0.0754	0.0705	0.0654	0.0637
170919	0.1527		0.0733	0.0657	0.0632	0.0620
246456	0.1273		0.0701	0.0632	0.0601	0.0589
554960	0.0849		0.0649	0.0553	0.0551	0.0550
1204362	0.0768		0.0610	0.0552	0.0550	0.0547

Equations

$${\Psi }_{\mathrm{pas}}(\mathbf{F})$$

$$\displaystyle\qquad+\dfrac{a_{fs}}{2b_{fs}}\bigl{[}e^{b_{fs}(I_{8,fs})^{2}}-1\bigr{]},$$

$${\mathbf{P}}_{\mathrm{pas}}$$

$$\mathbf{P}={\mathbf{P}}_{\mathrm{pas}}+{\mathbf{P}}_{\mathrm{act}},$$

$${\mathbf{P}}_{\text{act}}$$

$${\mathbf{\sigma }}_{\text{act}}$$

$$+\frac{{\kappa }_{nn}{T}_a}{J{\lambda }_n}\mathbf{F}{\mathbf{n}}_0\otimes \mathbf{F}{\mathbf{n}}_0,$$

$${\mathbf{F}}_A$$

$$+{\gamma }_s{\mathbf{s}}_0(\mathbf{x})\otimes {\mathbf{s}}_0(\mathbf{x})+{\gamma }_n{\mathbf{n}}_0(\mathbf{x})\otimes {\mathbf{n}}_0(\mathbf{x}).$$

$$\hat{\Psi }({\mathbf{F}}_E)$$

$$\displaystyle\quad+\dfrac{a_{fs}}{2b_{fs}}\bigl{[}e^{b_{fs}(I_{8,fs}^{E})^{2}}-1\bigr{]},$$

$$\mathbf{P}$$

$$I_1^E$$

$$-{\gamma }_n\frac{{\gamma }_n+2}{(1+{\gamma }_n{)}^2}{I}_{4,n},$$

$$I_{4,f}^E$$

$$I_{8,fs}^E$$

$$\det {\mathbf{F}}_A$$

$${\mathbf{F}}_A^{-1}$$

$$-\frac{{\gamma }_s}{1+{\gamma }_s}{\mathbf{s}}_0\otimes {\mathbf{s}}_0-\frac{{\gamma }_n}{1+{\gamma }_n}{\mathbf{n}}_0\otimes {\mathbf{n}}_0.$$

$$\dot{\mathbf{B}}=\nabla \mathbf{w}\mathbf{B}+\mathbf{B}(\nabla \mathbf{w}{)}^{\mathtt{t}},$$

$${\mathbf{\sigma }}_{\text{visc}}=\delta e^{\beta \dot{I_1}}\dot{\mathbf{B}},$$

$${\mathbf{P}}_{\text{tot}}=\mathbf{P}+J{\mathbf{\sigma }}_{\text{visc}}{\mathbf{F}}^{-\mathtt{t}},$$

$$\rho {\partial }_{tt}\mathbf{u}-\nabla \cdot {\mathbf{P}}_{\text{tot}}$$

$$\rho J-{\rho }_0$$

$$\mathcal{G}=\{\begin{array}{ll}\mathcal{G}(\mathbf{u},T_a)& :=\frac{\partial \Psi }{\partial \mathbf{F}}{\mathbf{F}}^{\mathtt{t}}+J{\mathbf{\sigma }}_{\text{visc}}+{\mathbf{P}}_{\text{act}}{\mathbf{F}}^{\mathtt{t}}\\ & \text{for active stress},\\ \mathcal{G}(\mathbf{u},\gamma )& :=\frac{\partial \hat{\Psi }}{\partial \mathbf{F}}{\mathbf{F}}^{\mathtt{t}}+J{\mathbf{\sigma }}_{\text{visc}}\\ & \text{for active strain},\end{array}$$

$$\mathbf{\Pi }=\mathcal{G}-pJ\mathbf{I}.$$

$$\mathbf{u}\cdot \mathbf{n}$$

$$\mathbf{\Pi }{\mathbf{F}}^{-\mathtt{t}}\mathbf{n}-p_{N}J{\mathbf{F}}^{-\mathtt{t}}\mathbf{n}$$

$$\mathbf{\Pi }{\mathbf{F}}^{-\mathtt{t}}\mathbf{n}+\eta J{\mathbf{F}}^{-\mathtt{t}}\mathbf{u}$$

$$\chi {\partial }_{t}v-\nabla \cdot \{\mathbf{D}(v,\mathbf{F},\mathbf{\Pi })\nabla v\}$$

$$\text{in }\Omega \times (0,t_{\mathrm{final}}],$$

$$\frac{d\vec{r}}{dt}$$

$$\text{in }\Omega \times (0,t_{\mathrm{final}}].$$

$$\mathbf{D}(v,\mathbf{F},\mathbf{\Pi })\nabla v\cdot \mathbf{n}$$

$$v(0)=0,\vec{r}$$

$$\mathbf{D}(v,\mathbf{F},\mathbf{\Pi })=$$

$$+D_{2}J{\mathbf{F}}^{-1}\mathbf{\Pi }{\mathbf{F}}^{-\mathtt{t}},$$

$${\partial }_t{T}_a=\hat{\alpha }\Delta T_a+\ell (T_a,\vec{r},I_{4,f})\text{in }\Omega \times (0,t_{\mathrm{final}}],$$

$${\gamma }_f(\xi )$$

$${\gamma }_n(\xi )$$

$$\frac{d\xi }{dt}=\hat{\ell }(\xi ,\vec{r})\text{in }\Omega \times (0,t_{\mathrm{final}}],$$

Full text

††journal: Biomech Model Mechanobiol

1 Introduction

In order to effectively combat cardiovascular disease, we need a robust scientific understanding of the mechanisms of the heart and the nature of such health conditions. Recent progress in the field is encouraging; the concept of patient-specific treatment is no longer a distant dream, but a conceivable reality and a topic of ongoing research. However, a major difficulty is our incomplete knowledge about the relationship between processes at the cellular and sub-cellular level, and the performance of the organ as a whole augustin16 . Indeed, a great deal of treatment is still based on trial-and-error experimentation rather than a more fundamental scientific understanding of the changes responsible for the onset and progression of disease goktepe13 . Several treatments, such as resynchronisation therapy and anti-arrhythmic medications, for example, are known to be ineffective or even exacerbate pathological conditions in some patients, for reasons that are not yet well understood jaffe14 . One potential obstacle to deep understanding of cardiac function is the difficulty of acquiring sufficiently detailed data. Until recently, there were no experimental techniques capable of recording 3D cardiac activity with high enough spatio-temporal resolution to provide the required level of information. However, relatively recent studies (see e.g. christoph17 ) have been able to use optical mapping to assess electromechanical waves with acceptable physiological accuracy.

Computational models have thus been critical in allowing for extensive study of the heart even without sufficient data. The development of complex multi-scale and multi-physics models, accompanied by advances in simulation and imaging techniques, has enabled researchers to investigate the many different aspects of cardiac function and disease. The hope is that the knowledge gained from these models can contribute to new and improved treatment methods. Even though the problem of cardiac electromechanics has been the focus of a large number of modelling and computational studies (see for instance augustin16 ; cherubini12 ; colli16 ; costabal17 ; gizzi15 ; goktepe13 ; nobile12 ; quarteroni17 ; sundnes14 and the references therein), there still remain many challenges in the development of more accurate and detailed models and the accompanying methods.

In such a context, the large majority of the proposed approaches rely on continuum formulations of the complex microstructural interactions occurring among the heart tissue components, e.g. cardiomyocytes, involving different scales quarteroni17 . The study of single cell and cell-cell lenarda18 chemomechanical and electromechanical interactions has attempted to unveil some of the underlying complex features of the cardiac function, and different multi-field nonlinear models have been gradually generalising classical approaches as the monodomain equations and Fick’s law of diffusion. In particular, fractional diffusion cusimano18 , nonlinear diffusion hurtado16 , and stress-assisted diffusion formulations cherubini17 were recently proposed to reproduce porous multiscale excitation phenomena within the framework of homogenised models for cardiac tissue. These studies, in fact, paved the route towards new challenging theoretical and computational problems aiming at a reliable in silico prediction of heart rate variability, cardiac repolarisation and inducibility of life-threatening arrhythmias phadumdeo18 . At the same time, macro-scale incompressibility, orthotropic and hysteretic mechanical features have been shown to fully characterise the human cardiac tissue under multiaxial loading tests (see e.g. gultekin16 and references therein). Viscosity properties, in particular, have been incorporated as one spring element coupled with Maxwell elements in parallel endowing the model with hysteretic characters describing the viscous response due to matrix, fibre, sheet and fibre-sheet couplings through four dedicated dashpots gultekin16 . Also in this case, a porous medium motivation has been advanced in yao12 , including the extracellular fluid filtrating through the elastic body, contributed by the active contractile behaviour of the muscle. However, complete agreement concerning the specific multiscale features involved in energy dissipation for the cardiac tissue, and soft biological tissues in general, is still lacking. The stress evolution equations for time-dependent viscous behaviour are based on finite strain viscoelasticity holzapfel01 , motivated by a rheological analogue from simo87 , and endowed with equilibrium and non-equilibrium contributions lubliner85 in which the usual assumption of volume-preserving deformations during time-dependent responses is made.

A distinguishing feature of our approach is the introduction of the mechanoelectrical feedback (MEF) in the electric conductivities, through a direct dependence on the Kirchhoff stress. This framework, known as stress-assisted diffusion (SAD), is widely employed in the modelling of gels and polymers klepach14 , but has only recently been adapted for active biological media undergoing reaction-diffusion excitation cherubini17 , and more tailored for cardiac models in loppini18 . While these contributions consider hyperelastic formulations coupled with multiphysics activation mechanisms, we also consider here the viscoelastic effects typical of soft microstructured fibre-reinforced biological tissues, and using realistic ventricular geometries. Fully mixed methods for the hyperelasticity of the cardiac tissue (that is, formulations involving stresses or strains in addition to simply displacement or displacement-pressure) are not yet widely employed. They have been introduced in ruiz15 and used more recently in cherubini17 ; garcia19 ; ruiz18 . In our case, our model and our numerical method include a three-field elasticity formulation (variationally based on a modification of the Hu-Washizu principle lami06 ) that states the governing equations in terms of stress-displacement-pressure, and that is motivated by the desire to avert volumetric locking and to solve directly for additional variables of interest. In particular, we solve for the Kirchhoff stress, which we use explicitly in our incorporation of SAD. This formulation includes a pressure-stabilisation term needed, in the lowest-order case, for triangular or for tetrahedral meshes. It constitutes a generalisation of the three-field formulation for nearly incompressible hyperelasticity, designed in chavan07 using quadrilateral meshes. Another difference in the present contribution is that we employ a more accurate cellular model, tailored for recovering human action potential dynamics, restitution features under constant pacing as well as sustained fibrillation behaviours and spiral waves breakup bueno08 . While the active strain approach is adopted in many instances in the literature and is often favoured due to the practicality of measuring strains directly using imaging techniques rossi14 , the active stress approach is somewhat simpler and more naturally incorporated in already existing models for passive deformation giantesio18 . In this work, we will adopt both formulations, although we find that the active strain formulation better reproduces physiologically accurate deformation regimes in ventricular geometries. To the best of our knowledge, no previous attempts have been made incorporating both active stress and active strain within a generalised stress-assisted reaction-diffusion formalism and embedding orthotropy, incompressibility and viscoelasticity for human cardiac ventricular domains.

This paper has been structured in the following manner. Section 2 lays out the elements of the mathematical model that describes the electro-viscoelastic function of the heart, including the active contraction of the cardiac muscle and the representation of the mechanoelectric feedback using stress-assisted conductivity, as well as a contribution from geometric nonlinearities (or geometrical feedback). The passive hyperelastic response of the tissue is described by an orthotropic exponential model, whereas the ionic activity which causes active contraction is incorporated through orthotropic active stress (active strain will also be addressed). The specific structure of the governing equations (written in terms of stress, displacements, electric potential, activation generation, and ionic variables) suggests to cast the problem in a mixed-primal form, and to use a mixed-primal finite element method for its numerical approximation. This is precisely the method that we outline in Section 3, which also includes a description of the consistent linearisation and implementation details. Our computational results in 2D and 3D, along with numerical validation and pertinent discussions on the modelling considerations are then presented in Section 4. We close with a summary and some remarks on model limitations and ongoing extensions, collected in Section 5.

2 Mathematical model

2.1 Finite-strain cardiac mechanics

Let $\Omega\subset\mathbb{R}^{d}$, $d\in\{2,3\}$ denote a deformable body with a piecewise smooth boundary $\partial\Omega$, considered in its reference configuration, and let $\boldsymbol{n}$ denote the outward unit normal vector on $\partial\Omega$. The kinematical description of finite deformations regarded on a time interval $t\in(0,t_{\text{final}}]$ is made precise as follows. A material point in $\Omega$ is denoted by $\boldsymbol{x}$, whereas $\boldsymbol{u}(t):\Omega\to\mathbb{R}^{d}$ will denote the displacement field defining its new position in the deformed configuration. The tensor $\mathbf{F}:=\mathbf{I}+\nabla\boldsymbol{u}$ is the gradient (applied with respect to the fixed material coordinates) of the deformation map; its determinant, denoted by $J=\det\mathbf{F}$, measures the solid volume change during the deformation; and $\mathbf{C}=\mathbf{F}^{\tt t}\mathbf{F}$ and $\mathbf{B}=\mathbf{F}\mathbf{F}^{\tt t}$ are respectively the right and left Cauchy-Green deformation tensors on which all strain measures will be based (here the superscript $()^{\tt t}$ denotes the transpose operator). The first isotropic invariant ruling deviatoric effects is the scalar $I_{1}(\mathbf{C})=\rm tr\,\mathbf{C}$, and for generic unitary vectors ${\boldsymbol{f}_{0}},{\boldsymbol{s}_{0}}$, the scalars $I_{4,f}(\mathbf{C})={\boldsymbol{f}_{0}}\cdot(\mathbf{C}{\boldsymbol{f}_{0}})$, $I_{4,s}(\mathbf{C})={\boldsymbol{s}_{0}}\cdot(\mathbf{C}{\boldsymbol{s}_{0}})$, $I_{8,fs}(\mathbf{C})={\boldsymbol{f}_{0}}\cdot(\mathbf{C}{\boldsymbol{s}_{0}})$ are pseudo-invariants of $\mathbf{C}$ measuring direction-specific stretch ciarlet .

The triplet $({\boldsymbol{f}_{0}}(\boldsymbol{x}),{\boldsymbol{s}_{0}}(\boldsymbol{x}),{\boldsymbol{n}_{0}}(\boldsymbol{x}))$ represents a coordinate system pointing in the local direction of the muscular cardiac fibres, transversal sheetlet compound, and normal cross-fibre direction ${\boldsymbol{n}_{0}}(\boldsymbol{x})={\boldsymbol{f}_{0}}(\boldsymbol{x})\times{\boldsymbol{s}_{0}}(\boldsymbol{x})$. Note that the system is restricted to $({\boldsymbol{f}_{0}}(\boldsymbol{x}),{\boldsymbol{s}_{0}}(\boldsymbol{x}))$ in the two-dimensional case, and that these directions are defined in the reference configuration. Constitutive relations characterising the material properties and underlying microstructure of the myocardial tissue will follow the orthotropic model proposed in holzapfel09 , whose strain energy density (relating the amount of energy stored within the material in response [Joule/Volume] to strain, and which assumes an additive decomposition into isotropic, volumetric and anisotropic contributions) and the first Piola-Kirchhoff stress tensor (associated with a passive, elastic deformation) read, respectively

[TABLE]

where $a,b$ are material constants associated with the isotropic matrix response, $a_{f}$ and $b_{f}$ rule the directional behaviour of the material along myocardial fibres, $a_{s}$ and $b_{s}$ account for the cross-contribution of the fibre sheet directions, and $a_{fs},b_{fs}$ encapsulate the shear effects in the fibre-sheet plane. Moreover, the field $p$ denotes the solid hydrostatic pressure, and we have used the notation $(u)_{+}:=u$ if $u>0$ or zero otherwise, for a generic real-valued function $u$. This modelling choice aligns with the fact that fibres have a quite different behaviour under compression or tension regimes. In addition, taking the positive part of the exponents in the anisotropic energy results in excluding anisotropic energetic contributions for compressed fibre configurations, which in the case of passive fibres should have an effect only during extension pezzuto14 . We remark here that the particular mechanisms of soft tissue anisotropic mechanical behaviour is still under investigation humphrey14 . Moreover, full incompressibility of the tissue will be enforced in the present framework, and this has some advantages associated with the mathematical and numerical structure of the system. Although biological tissues possess a complex porous structure, compression features are still being systematically investigated ex-vivo, and a more comprehensive answer on the subject is still needed mcevoy18 .

2.2 Active stress and active strain

In physiological scenarios, the mechanical deformation is also actively influenced by microscopic tension generation.

Active stress model. A simple description is given in terms of active stresses (see for instance sundnes14 ): we assume that the first Piola-Kirchhoff stress tensor decomposes as

[TABLE]

where the active stress component acts differently on each local direction with an intensity depending on the scalar field of active tension $T_{a}$, that synthesises the biochemical state of myocytes (and whose dynamic behaviour will be specified later on). Then

[TABLE]

where $\kappa_{sn},\kappa_{nn}$ are positive constants representing the variation of activation on each specific direction, as proposed in dorri06 , and $\lambda_{f}=\sqrt{I_{4,f}},\lambda_{s}=\sqrt{I_{4,s}},\lambda_{n}=\sqrt{{\boldsymbol{n}_{0}}\cdot(\mathbf{C}{\boldsymbol{n}_{0}})}$ are the fibre, sheetlet, and cross - fibre stretches. Setting appropriate models for $\boldsymbol{\sigma}_{\text{act}}$ is not a trivial task since the active contribution to the force should account for the geometric properties of deformation, and these undergo substantial changes during contraction in the finite strain regime pezzuto14 . Details of other anisotropic activation forms can be found, for instance, in rossi14 for active strain and in (usyk00, , Appendix B) for active stress descriptions, but they are basically responsible for additional deformation effects such as wall thickening, radial constriction and torsion, as well as longitudinal shortening. Note that the active Cauchy stress does not include a contribution on the diagonal entry associated with the local sheetlet direction ${\boldsymbol{s}_{0}}$ since a stress component on this direction would counteract wall thickening mechanisms dorri06 . Moreover, the intensity of the active tension effect on the cross-fibre direction ${\boldsymbol{n}_{0}}$ is assumed to be substantially smaller than that appearing on the off-diagonal component $\mathrm{sym}(\mathbf{F}{\boldsymbol{s}_{0}}\otimes\mathbf{F}{\boldsymbol{n}_{0}})$, see Table 4.2. Also note that some references do not include a rescaling with local stretches in each term of $\boldsymbol{\sigma}_{\text{act}}$.

Active strain model. Next we recall the active strain model for ventricular electromechanics (see e.g. cherubini08 ). There, the contraction of the tissue results from activation mechanisms governed by internal variables and incorporated into the finite elasticity context using a multiplicative decomposition of the deformation gradient into a passive (purely elastic) and an active part, $\mathbf{F}=\mathbf{F}_{\!E}\mathbf{F}_{\!A}$, with

[TABLE]

The coefficients $\gamma_{i}$, with $i=f,s,n$, are smooth scalar functions encoding the macroscopic stretch in specific directions, whose precise definition will be postponed to (2.15). The inelastic contribution to the deformation modifies the length and maybe also the shape of the cardiac fibres, and then compatibility of the motion is restored through an elastic deformation accommodating the active strain distortion. A physiological motivation for the active strain approach is related to the shortening of sarcomeres as a response to the sliding filaments of the actin-myosin molecular motor: such shortening is encapsulated in $\mathbf{F}_{\!A}$, which determines a new (and fictitious, or virtual) intermediate configuration that is regarded as a reference for the elastic deformation pezzuto14 . Therefore, the strain energy function and the first Piola-Kirchhoff stress tensor (after applying the active strain decomposition) are functions of $\mathbf{F}_{\!E}$ only, and they read respectively

[TABLE]

As in the description of (2.1) above, we again note that one switches off the anisotropic contributions under compression. An additional advantage is that the associated terms in the strain energy function (in both the pure passive and active-strain formulations) can be shown to be strongly elliptic pezzuto14 (these will be the terms appearing on the second diagonal block of the weak formulation from Section 3, the block corresponding to displacements), however the overall problem will remain of a saddle-point structure. The modified elastic invariants $I_{i}^{E}$ are functions of the coefficients $\gamma_{i}$, as well as of the invariant and pseudo-invariants in the following manner rossi12

[TABLE]

Such dependencies are a consequence of assuming isochoric active deformations pezzuto14 , i.e. $\det\mathbf{F}_{\!A}=1$, justified by the fact that the volume of the cardiomyocytes does not vary substantially during contraction. Besides, following Rossi et al. rossi12 , previous expressions are obtained by assuming $\gamma_{s}=\gamma_{n}$ and making use of the fact that $I_{1}=I_{4,f}+I_{4,s},I_{4,n}$, as well as that $\mathbf{F}_{\!E}=\mathbf{F}\mathbf{F}_{\!A}^{-1}$, with

[TABLE]

Accordingly, the active strain, and consequently the force associated to the active part of the total stress, will receive contributions acting on the three main directions. The calcium-based activation signal travels up to four times faster along the fibre axis than in the sheet and normal directions, and this fact further motivates the use of orthotropic active strain rossi14 .

2.3 Viscoelasticity and equations of motion

Extension and shear tests demonstrate the importance of incorporating viscoelastic effects in models for cardiac passive mechanics gultekin16 . In the heart, the extracellular fluid filtrating through the elastic solid is one of the main generators of the viscoelastic effects of the tissue yao12 . Viscous effects are also tied to crossbridge processes identified in Ca2+ activated fibres maughan98 , and have a well-established literature as well as a consistent methodology for their implementation (the stress update algorithm that uses a convolution integral representation) developed for general soft tissues holzapfel01 . From the viewpoint of kinematics, it suffices to relate stress to strain rates. Decomposition of the spatial velocity gradient $\boldsymbol{w}=\dot{\boldsymbol{u}}$ into the rate of deformation and spin tensors yields the relation

[TABLE]

and a simplified rheological Kelvin-Voigt model for the viscous component of the Cauchy stress can be defined as follows (see e.g. karlsen17 )

[TABLE]

which depends on the history of the isotropic contribution to the Cauchy stress. Here $\delta,\beta>0$ are model parameters. In this way, after a pull-back operation, we see that

[TABLE]

is the total first Piola-Kirchhoff stress tensor that includes $\mathbf{P}$ defined from either (2.1)-(2.2) or (2.2), and the viscoelastic contributions.

More advanced rheologies can be easily incorporated in the context of active stress formulations as done in e.g. katsnelson04 , as the generalised Hill-Maxwell model recently proposed in cansiz17 , as in the perturbed equations of harmonic wave motion using springpot-based models with fractional order derivatives from capilnasiu19 , or as the thermodynamical electro-viscoelastic models that use statistical fibre distributions gizzi18b . We will, however, confine the presentation to (2.7) without introducing stochasticity of the anisotropic components.

Irrespective of the activation formalism one adopts (active strain or active stress), the balance of linear momentum and the incompressibility constraint (allowing only isochoric motions) are written together in the following way, when posed in the inertial reference frame and under transient mechanical equilibrium,

[TABLE]

where $\rho_{0},\rho$ are the reference and current medium density, $\boldsymbol{b}$ is a smooth vector field of imposed body loads, $\partial_{tt}$ denotes the second time derivative, and the divergence operator in (2.8a) applies on the tensor fields row-by-row. The balance of angular momentum translates into the condition that the Kirchhoff stress tensor $\boldsymbol{\Pi}=\mathbf{P}_{\text{tot}}\mathbf{F}^{\tt t}$ must be symmetric, which is in turn encapsulated into the momentum and constitutive relations (2.8a), (2.1), (2.2).

Defining

[TABLE]

as the contribution to the Kirchhoff stress that does not involve pressure, we then have

[TABLE]

Stating the balance equations in terms of Kirchhoff stress, displacements, and pressure suggests that, at the level of writing finite element schemes, we will use mixed methods respecting the same structure. Additionally, setting boundary conditions for the motion of the left ventricle is not trivial, as the organ is known to slightly move and twist during the heartbeat. In our case, equations (2.8a)-(2.8b)-(2.9) will be supplemented with mixed normal displacement and traction boundary conditions

[TABLE]

where $\partial\Omega_{D}$, $\partial\Omega_{N}$, $\partial\Omega_{R}$ conform a disjoint partition of the boundary. The condition (2.10a) means that we constrain the normal motion along the normal direction with respect to the surface $\partial\Omega_{D}$. The term $p_{N}$ in (2.10b) denotes a (possibly time dependent) prescribed boundary pressure associated with endocardial blood pressure, which is uniform over the deformed counterpart of $\partial\Omega_{N}$ and it is applied in the normal direction to the epicardium in the deformed configuration. However, owing to Nanson’s formula ciarlet , this contribution regarded on the reference configuration depends on the cofactor of the deformation gradient and therefore the boundary condition is nonlinear in the undeformed configuration; moreover, the traction written in terms of the Kirchhoff stress tensor is $\boldsymbol{t}=\boldsymbol{\Pi}\mathbf{F}^{-\tt t}\boldsymbol{n}$. Also note that the Robin conditions (2.10c) account for stiff springs connecting the cardiac medium with the surrounding soft tissue and organs (whose stiffness is encoded in the scalar $\eta$). More sophisticated boundary conditions that consider an interaction with the pericardium can be also imposed fritz14 . A sketch of a mono-ventricular domain specifying boundary surfaces and fibre directions is depicted in Figure 2.1.

2.4 Monodomain equations

In the context of electromechanical processes, the propagation of electric potential $v$ is governed by the following reaction-diffusion system, known as the monodomain equations (see e.g. ColliBook ), which are cast here in the reference configuration. The current conservation is written only in terms of the transmembrane potential and the coupling with additional ionic quantities are encoded in the vector $\vec{r}$ (here we use $\vec{\cdot}$ instead of bold to denote vector fields of dimension other than $d$)

[TABLE]

Here $\chi$ is the ratio of membrane area per tissue volume, and $I_{\rm ext}$ is a spatio-temporal external stimulus applied to the medium. We will adopt the minimal model for human ventricular action potential, proposed in bueno08 and fitted to capture restitution curves, conduction velocity, spiral/arrhythmic dynamics, and complex behaviour typical to nonlinear dynamical systems (used later for cardiac alternans in gizzi13 ). That model was, however, tailored originally for the case of isotropic conductivity $\mathbf{D}=D\mathbf{I}$, and so the extended fully-coupled model discussed below will be able to accommodate a wider class of propagation patterns, and will also constitute a generalisation over other recent models for stress-assisted diffusion cherubini17 ; loppini18 .

Specification of the ionic currents and gating variables can be found in Appendix A.

Boundary and initial conditions for (2.11) correspond to

[TABLE]

and (2.12b) can be combined with suitable initial pacing, especially needed in more complex and more physiologically accurate cell models. The minimal model, as proposed in bueno08 , has a heterogeneous character that we do not consider in our study. Their description contains separate parameter sets that are able to reproduce experimental results for the epicardium, mid-myocardium and endocardium, as well as parameter sets that mimic the results of two more complicated ionic models for human ventricular cells. For simplicity (and also as a consequence of lack of personalised experimental data) we use the parameter set developed for the epicardium (see values in Table 4.2), assuming that it is consistent throughout the cardiac wall. Extension to the heterogeneous case can be readily incorporated.

2.5 Stress-assisted conduction

The mechanoelectrical feedback (the process where the current mechanical state of the deforming solid modifies both the excitability and electrical conduction of the tissue) is here introduced in the conductivity tensor, through a direct dependence on the Kirchhoff stress (which constitutes one of the main novelties in our approach, stemming as a generalisation of the anisotropy induced by stress proposed in cherubini17 and later used for simplified 2D cardiac electromechanics in loppini18 ). In addition, due to the Piola transformation (yielding a transformation of the diffusion tensor using the deformation gradients), the conductivity tensor also depends nonlinearly on the deformation gradient (actually, the term $J\mathbf{C}^{-1}$ constitutes a strain-based modification of tissue conductivity, also referred to as geometric feedback in colli16 )

[TABLE]

where the nonlinear conductivity (self diffusion depending on $v$) accounts for porous media electrophysiology following the development in hurtado16 , but appropriately modified to incorporate information about preferred directions of diffusivity according to the microstructure of the tissue (encoded in the second term defining $\mathbf{D}$). The parameter $D_{0}$ signifies the usual diffusion for isotropic materials, whereas $D_{1}$ and $D_{2}$ represent the intensity of the porous media electrophysiology and of the stress-assisted diffusion, respectively. An additional term in the nonlinear self-diffusion (e.g. $D_{3}v^{2}$, as in gizzi17 ; ruiz18 ) eventually leads to very slight modifications in conduction patterns and we have therefore decided not to include it. Tuning $D_{1}$ is sufficient to, if needed, calibrate the speed and action potential duration at the depolarisation plateau phase.

It is useful to point out that both the nonlinear self-diffusion term and the SAD argument derive from rigorous thermodynamical principles, formulated under specific assumptions for porous materials. In particular, nonlinear self-diffusion is naturally related to the transport of chemicals within porous media, while classical models of stress-assisted diffusion for general materials aifantis80 also consider the transport of diffusing species within solids exhibiting finite strains. For the specific case of cardiac tissue, both approaches are justified by the multiple scales involved in the transport of ions and generation and propagation of action potential within the cell and across different cells lenarda18 . In particular, we can mention the role of intercalated discs and gap junctions between communicating cells or the presence of the ephathic couplings in the extracellular space ly18 , as well as micro-invaginations on the cell membrane known as microtubules and microdomains miragoli16 . All of these emerging effects contribute to the macroscopic nonlinearities and additional anisotropies considered in the diffusion tensor herein and which could be further analysed through a consistent multiscale homogenisation study.

It is important to remark that the solvability of the monodomain equations (2.11a)-(2.11b) depends on the properties of $\mathbf{D}$. In particular, the stress-assisted diffusion tensor needs to remain symmetric and uniformly elliptic, which is a non-trivial condition, given the dependence on stress and on voltage. A thorough sensitivity analysis (but for a simpler dependence on stress) can be found in cherubini17 . Here we perform a much lighter calibration, as mentioned later in Section 4. Comparisons between the effects of SAD and the more conventional mechanoelectrical feedback through stretch-activated currents have been reported in loppini18 .

2.6 Activation and excitation-contraction coupling

When using the active stress approach, we will adopt a simple description where the active tension is generated by ionic quantities (calcium) as well as by local fibre stretch. That is, we propose a regularised active tension model of the form

[TABLE]

with $\hat{\alpha}=\alpha_{1}D_{0}$, and $\ell(T_{a},\vec{r},I_{4,f})=T_{a}-\alpha_{2}r_{3}+\alpha_{3}I_{4,f}$, where $\alpha_{1},\alpha_{2},\alpha_{3}=0.1\alpha_{2}$ are positive model constants. As calcium concentration is not readily available in the phenomenological cellular model we are employing, we use $r_{3}$ as a proxy for intracellular calcium bueno08 . In addition, a linear dependence on the calcium proxy and on the local stretch are sufficient in our setting to qualitatively capture the dynamics of active tension.

On the other hand, in the framework of active strain, a constitutive equation for the activation functions $\gamma_{i}$ in terms of the microscopic cell shortening $\xi$ is expressed as follows (see e.g. barbarotta18 )

[TABLE]

and the specific relation between the myocyte shortening $\xi$ and the dynamics of slow ionic quantities (in the context of our phenomenological model, only $\vec{r}$) is made precise using the law

[TABLE]

which does not require an explicit dependence on local fibre stretch, as the sliding of myofilaments is driving the dynamics of the functions $\gamma_{i}$. We employ the nonlinear reaction term $\hat{\ell}(\xi,\vec{r})=K_{1}(1+r_{3})^{-1}+K_{2}\xi$, and we make the distinction that $\ell$ and $\hat{\ell}$ characterise the evolution of the activation in the approaches of active stress and active strain, respectively.

3 Numerical method and implementation

3.1 Mixed-primal weak form

Restricting to the case of an active strain model with Robin conditions (2.10c) on the whole boundary for the mechanical layer (that is $\partial\Omega_{R}=\partial\Omega$) and the boundary and initial conditions (2.12a)-(2.12b) for the electrical layer, we proceed to take the inner product of the differential equations (2.8a), (2.8b), (2.9), (2.11), (2.16) with adequate test functions, and to integrate by parts whenever appropriate. We then arrive at the following weak form of the problem: For $t>0$, find $(\boldsymbol{\Pi},\boldsymbol{u},p)\in\mathbb{L}^{2}_{\text{sym}}(\Omega)\times\mathbf{H}^{1}(\Omega)\times\mathrm{L}^{2}(\Omega)$ and $(v,\vec{r},\xi)\in\mathrm{H}^{1}(\Omega)\times\mathrm{L}^{2}(\Omega)^{3}\times\mathrm{L}^{2}(\Omega)$ such that

where $\mathbb{L}^{2}_{\text{sym}}(\Omega):=\{\boldsymbol{\tau}\in\mathbb{L}^{2}(\Omega):\boldsymbol{\tau}=\boldsymbol{\tau}^{\tt t}\}$, and where the case for an active stress formulation necessitating an active tension model is addressed similarly (however, the regularity of $T_{a}(t)$ is then $\mathrm{H}^{1}(\Omega)$). Theoretical aspects regarding the coupling of elasticity and stress-assisted diffusion problems has been recently addressed in the context of mixed-primal and mixed-mixed formulations in gatica18b , but only for the case of simplified linear three-field elasticity and steady diffusion.

The spatial discretisation will follow a mixed-primal Galerkin approach based on the weak formulation (LABEL:eq:weak-form). Details on the finite-dimensional spaces and linearisation procedure are laid out in Appendix B.

The motivation for using three-field elasticity formulations is the need to produce robust solutions with balanced convergence orders for all variables. In addition, these methods are robust in the incompressible regime; they are not subject to volumetric locking lamichhane12 ; and most importantly, they provide direct approximation of variables of interest, albeit at a higher computational cost. Another advantage of using the Kirchhoff stress is that this tensor is symmetric, and, for simpler material laws, is a polynomial function of the displacements (whereas first and second Piola-Kirchhoff stresses are rational functions of displacement) chavan07 . Solving in terms of stresses proves particularly useful, as this variable participates actively in the electromechanical coupling through the stress-assisted diffusion. Moreover, for the lowest-order method characterised by $l=0$, the matrix system associated with (LABEL:eq:FE) has fewer unknowns than the discretisation that uses piecewise quadratic and continuous displacement approximations and piecewise linear and discontinuous pressure approximations (and which is a popular locking-free scheme for hyperelasticity in the displacement-pressure formulation, utilised for stress-assisted diffusion problems in the recent work loppini18 ). The importance of casting the equations of motion in terms of the coupling variables has been already emphasized in ruiz15 in the context of cardiac electromechanics, which demonstrates that the computation of output indicators of interest (such as conduction velocities) may suffer from loss of accuracy if one simply postprocesses stress or strain from discrete displacements as approximations in the geometric feedback.

3.2 Solver structure and implementation details

According to the fixed-point separation between electrophysiology and solid deformation solvers, nonlinear mechanics will be solved using the Newton-Raphson method stated above, and an operator splitting algorithm will separate an implicit diffusion solution (where another Newton iteration handles the nonlinear self-diffusion) from an explicit reaction step for the kinetic equations, turning the overall solver into a semi-implicit method. Such a strategy is feasible since the Jacobians associated with the reaction and excitation-contraction models do not possess highly varying eigenvalues (otherwise one would need to include these terms in the Newton iteration). Updating and storing of the internal variables $\xi$ and $\vec{r}$ will be done locally at the quadrature points. We solve the linear systems arising at each Newton iterate by the Krylov iterative method GMRES, preconditioned with an incomplete LU(0) factorisation (except for the linear systems in the convergence tests in Section 4.1, which will be solved with the direct method SuperLU), and the iterates are terminated once a tolerance of $10^{-6}$ (imposed on the $\ell^{\infty}$ norm of the non-preconditioned residual) has been achieved. The mass matrices associated with the discretisation of the monodomain equations are assembled in a lumped manner, which reduces the amount of artificial diffusion and violation of the discrete maximum principle quarteroni17 . All routines have been implemented using the finite element library FEniCS fenics .

4 Computational results

4.1 Mesh convergence

We begin with the numerical validation of our mixed-primal method on a problem slightly simpler than (2.8) - (2.11) - (2.14), but that still retains the main ingredients of the model. These include orthotropic active mechanics, nonlinear reaction-diffusion with stress-assisted diffusion, and a nonlinear excitation-contraction coupling.

A convergence test is generated by computing errors between smooth exact solutions and approximate solutions using the first-order and the second-order methods discussed in Section 3. Let us consider the following closed-form solutions to a steady-state counterpart of the variational form (LABEL:eq:weak-form) for the electromechanics equations, also assuming the absence of viscoelastic effects, and defined on the domain $\Omega=(0,1)^{2}$ with the fibres/sheetlets defined as ${\boldsymbol{f}_{0}}=(0,1)^{\tt t},{\boldsymbol{s}_{0}}=(-1,0)^{\tt t}$

[TABLE]

Then the Kirchhoff stress $\boldsymbol{\Pi}$, as well as suitable forcing terms (volumetric load, an additional external stimulus, and the active tension source) are computed from these smooth solutions, the balance equations, relations (2.2), (2.3), (2.13), and using the following simplified constitutive equations

[TABLE]

Note also that the incompressibility constraint for this test is $J=J_{ex}$, where $J_{ex}$ is computed from the exact displacement. Here we also prescribe Dirichlet boundary conditions for displacements, transmembrane potential, and active tension (incorporated in the discrete trial spaces). Errors due to fixed-point iterations are avoided by taking a full monolithic coupling and computing solutions using Newton-Raphson iterations with an exact Jacobian. On a sequence of six uniformly refined meshes, we proceed to compute errors between the exact and approximate solutions computed with methods using $l=0$ and $l=1$. Kirchhoff stress and pressure errors are measured in the $L^{2}-$norm, whereas for the remaining variables the errors are measured in the $H^{1}-$norm. The obtained error history is reported in Table 4.1, where we observe an asymptotic $O(h^{l+1})$ decay of the error for each field variable. This behaviour corresponds to the optimal convergence according to the interpolation properties of the employed finite element subspaces chavan07 .

4.2 Parameter calibration

For the following 2D simulations, we will initially consider tissue slabs of $50\times 50$ mm2, and set fibre and sheetlet directions simply as ${\boldsymbol{f}_{0}}=(1,0)^{\tt t}$, ${\boldsymbol{s}_{0}}=(0,-1)^{\tt t}$. The initiation, maintenance, prevention and treatment of so-called reentrant waves is a major focus of current research due to their implication in atrial and ventricular fibrillations ColliBook . We are thus interested in investigating the formation of spiral reentrant waves in our model setup, following the S1-S2 stimulation protocol. We first excite the tissue with a symmetric stimulus labelled S1. An asymmetric stimulus labelled S2 is then applied during the vulnerable window near the end of the refractory period, when some of the tissue has recovered excitability but depolarisation is still blocked elsewhere. This causes unbalanced excitation, which can lead to the formation of a spiral wave. We will define the spiral front as the edge of the spiral wave, where the excitation front meets the repolarisation waveback of the action potential. In our simulations, both waves have nondimensional amplitude 3 and duration 3 ms. The S1 stimulus is a planar wave created by exciting the entire left edge of the tissue in 2D, and the entire bottom section (below some value of the $z-$coordinate) in 3D. The S2 stimulus is a square wave created by exciting the bottom left quadrant at $t=330$ ms and the bottom left octant at $t=335$ ms in 2D and 3D, respectively. Here we use the active stress approach, and the boundary conditions for the visco-elastodynamic equations correspond to (2.10c). The formation and evolution of the spiral wave on a deforming domain can be seen in Figure 4.1. The spiral is initiated by the diffusion of voltage and transport of ionic entities from the S2 stimulus into the leftmost section of the tissue, which has recovered enough excitability after S1. The wave then spreads outwards in all directions, occupying the entire tissue except for the region that was just excited by the S2 wave.

Next, since we are using the active stress formulation in this case, we proceed to evaluate $\alpha_{1},\alpha_{2}$, the parameters governing active tension in (2.14), and $\eta$, the stiffness parameter from (2.10c). We conduct a simple sensitivity analysis by increasing or decreasing either $\alpha_{1},\alpha_{2}$ or $\eta$ by one order of magnitude, holding the others constant at their reference values ($\alpha_{1}=10$, $\alpha_{2}=0.5$, and $\eta=\eta_{a}=0.001$N$/$cm2, as listed in Table 4.2). This simple analysis therefore does not test for compounding or interaction effects. We also consider a smaller slab of size $12\times 12$ mm2. The parameter $\alpha_{1}$ contributes to producing smoother active tension profiles, while $\alpha_{2}$ controls the range of their magnitude. These effects are visible in Figure 4.2. We found that larger values of $\alpha_{1}$ produced smoother gradients in pressure and stress, while larger values of $\alpha_{2}$ produced, in average, higher magnitude displacement, Kirchhoff stress, and pressure, as well as some more subtle changes in ionic quantities. Parameter $\eta$ determines the stiffness of the springs supporting the tissue, and so decreasing $\eta$ resulted in an increase in the maximum values of magnitude of displacement, stress, and pressure, as expected. However, these differences were minimal, even across the three orders of magnitude tested ($\eta=1$E-4 to $\eta=0.01$). The effects on ionic entities were even smaller, for both the hyperelastic and viscoelastic cases, and therefore plots are not shown.

Computational experiments reveal a window of values of $D_{2}$ for which our method converges. In the 2D hyperelastic case, we found that the upper bound for $D_{2}$ is approximately $D_{2}=2.1$E-2 mm${}^{2}/$s, with the linear solver failing to converge for larger values. In these simulations, the Kirchhoff stress achieved an $L^{2}-$norm of between 0.006 and 0.6. In turn, the viscoelastic case was able to accept slightly larger values of $D_{2}$, up to $D_{2}=2.2$E-2 mm${}^{2}/$s, with the $L^{2}-$norm of stress falling between 0.001 and 0.5. A possible explanation is the loss of coercivity or monotonicity in the stress-assisted diffusion coupling, as explored in cherubini17 .

The numerical method used for these tests is characterised by the time step, mesh size, polynomial degree, and stabilisation constant $\Delta t=0.01$ ms, $h=0.3534$ mm, $\ell=0$, $\zeta_{\mathrm{stab}}=2.5$, respectively.

4.3 Locking-free property

We next proceed to assess the performance of the proposed mixed formulation for the mechanical problem. In this example we solve only for (2.8) without the acceleration term (otherwise present in all other simulations), using the active stress approach with a fixed value for the active tension and without the contribution from the viscous stress (2.6). We proceed to compare the deformation achieved by the mixed formulation with that of an asymptotic solution and the approximate solution generated by a more standard pressure-displacement finite element formulation. We consider different stabilisation parameter values and mesh refinements.

We perform two sets of computations. First, we undertake Cook’s membrane benchmark test for a fully incompressible Holzapfel-Ogden material (as was similarly done for nearly incompressible Saint Venant-Kirchhoff and Neo-Hookean solids in (chavan07, , Test II)), where we set an active tension of $T_{a}=0.07$. This test involves applying an upward in-plane shear load $\boldsymbol{t}=(0,100)^{\tt t}$ to the right edge of a tapered panel with a clamped left edge, and measuring the vertical deformation of the upper right vertex. The domain is defined as the convex hull of the set $\{(0,0),(48,44),(48,60),(0,44)\}$ (see the sketch in Figure 4.3(a)), and the fibre and sheetlet fields are ${\boldsymbol{f}_{0}}=(1,0)^{\tt t}$ and ${\boldsymbol{s}_{0}}=(0,-1)^{\tt t}$, respectively. Secondly, we consider a 3D system suggested in (land15, , Test I) as a simple benchmark for passive cardiac mechanics, and therefore we set $T_{a}=0$. The problem consists in computing the deformation of a point at the right end of a beam defined by the domain $\Omega=(0,10)\times(0,1)\times(0,1)$  mm (see the sketch in Figure 4.3(b)), where the fibre direction is ${\boldsymbol{f}_{0}}=(1,0,0)^{\tt t}$. Instead of (2.1), the material is characterised by the transversally isotropic strain energy density proposed by Guccione et al. gucc (which is the material law used in the benchmark test from land15 ): $\Psi_{\text{pas}}=a/2(e^{Q}-1)$, with $Q=b_{f}E_{ff}^{2}+b_{t}(E_{ss}^{2}+E_{nn}^{2}+E_{sn}^{2}+E_{ns}^{2})+b_{fs}(E_{fs}^{2}+E_{sf}^{2}+E_{fn}^{2}+E_{nf}^{2})$, where $a=2$ kPa, $b_{f}=8$, $b_{t}=2$, $b_{fs}=4$, and the $E_{ij}$ denote entries of the Green-Lagrange strain tensor $\mathbf{E}$, rotated with respect to a local coordinate system aligned with ${\boldsymbol{f}_{0}},{\boldsymbol{s}_{0}},{\boldsymbol{n}_{0}}$. The beam is clamped at the face $x=0$, a pressure of $p_{N}=0.004$ kPa is imposed on the bottom face $z=0$, and the remainder of the boundary is considered with traction-free conditions. According to (2.10b), the pressure boundary condition changes with the deformed surface orientation, and its magnitude scales with the deformed area.

The outcome of these tests in Figures 4.4(a,b) shows a rapid convergence of our first- and second-order methods, while the computations using a pressure-displacement formulation and the Taylor-Hood finite elements (the well-known $\mathbb{P}_{2}-\mathbb{P}_{1}$ pair of continuous and piecewise quadratic approximations of displacements and continuous and piecewise linear approximations for pressure) display a somewhat slower convergence to the asymptotic deflection of the membrane. Using discontinuous pressures (the $\mathbb{P}_{2}-\mathbb{P}_{1}^{\text{disc}}$ pair) rectifies the convergence, but at a higher computational cost. Quite similar results were obtained for the beam (where the reference value is the average of the reported simulations from the study in land15 ). Moreover, Figures 4.4(c,d) show the vertical deflections as a function of the number of vertices discretising the left side of the membrane and of the small edge of the beam, respectively. They indicate that the obtained results are consistent for varying values of the stabilisation parameter, $\zeta_{\text{stab}}$, and the observed behaviour also confirms that our method is locking-free.

4.4 Stress-assisted diffusion and conduction velocity

In addition to determining a suitable parameter range for $D_{2}$ that ensures solvability of the discrete monodomain equations, we also investigated the effect of $D_{2}$ on the tissue’s response to spiral wave dynamics. As in the second part of Section 4.2, this time the domain is a square slab of width $12\,{\rm mm}$ aligned with the canonical axes. We employ the active stress approach and use Robin boundary conditions for the viscoelasticity problem. The fibres assume the fixed direction ${\boldsymbol{f}_{0}}=(1,0)^{\tt t}$ and the sheetlets ${\boldsymbol{s}_{0}}=(0,-1)^{\tt t}$, and the Holzapfel-Ogden material law is considered. The meshsize is approximately $h=0.085\,$mm and the timestep is $\Delta t=0.1\,$ms. We use the lowest-order finite element method and the stabilisation parameter is $\zeta_{\text{stab}}=2.5$.

Figure 4.5 shows the differences in the ionic quantities between simulations with a very small contribution of SAD ($D_{2}=$ 1 E-5 mm${}^{2}/$s) and a more prominent, but still mild SAD contribution ($D_{2}=$7.5E-3 mm${}^{2}/$s). The snapshots correspond to the time $t=444$ ms, when the spiral tip has not yet formed. A closer inspection suggests that these contrasts were due to a difference in conduction velocity induced by SAD. In Figure 4.6(a,b), we see that conduction velocity was higher for larger values of $D_{2}$ (meaning a larger SAD contribution). When the wave first emerged, the peak action potential was more advanced for the case of reduced $D_{2}$, but the large $D_{2}$ peak eventually caught up to and surpassed it, which is a phenomenon also observed in the active tension curves. The ionic quantities followed the same trend. Indeed, an analysis similar to that which produced Figure 4.6(a,b) revealed that the overall profiles of the ionic quantities were highly similar between the two cases compared in Figure 4.5, but differed in the speed at which they are transported through the tissue.

We also remark that the effect of changing conduction velocities was not spatially consistent. SAD increases conduction velocity in the fibre (horizontal) direction, but actually decreased conduction velocity in the vertical and diagonal directions. This resulted in a noteworthy effect on the growth of the spiral wave. Figures 4.6(c,d,e) show a comparison of the spiral wave in the viscoelastic case for three different values of $D_{2}$. The upper right area of the spiral is slightly more vertical in the simulation with a larger value of $D_{2}$ than in the other cases, suggesting that propagation of the voltage is somewhat hindered in the fibre direction. We also observe a slightly more pronounced deformation of the right side of the domain due to the two-way coupling between tissue motion and electrophysiology. A similar effect was seen in the hyperelastic case.

As in other studies, here we observe that conduction velocity is sensitive to spatio-temporal discretisation. In Table 4.3, we include the results of a simple convergence test for conduction velocity, similar to the benchmark test conducted in ruiz18 . We calculated the horizontal propagation of the action potential using different time steps and mesh refinements. Differently to the case of nonlinear diffusion without SAD from ruiz18 , the experiment reveals that lower resolutions produce larger conduction velocities than the physiological values. This test also confirms that with our time step and mesh resolution ($0.1\,\text{ms}$, and above $200,000$ DoF, respectively), conduction velocity is in the expected physiological range, whereas larger time steps will systematically fail to capture the dynamics of the ionic model.

4.5 Scroll waves on mono-ventricular geometries

For the ventricular geometries, we test both the active strain and active stress formulations. We start from patient-specific left ventricular geometries (available from warriner18 ; lamata18 ) and rescale them using approximately the same dimensions as the idealised ventricles studied in ruiz18 . The segmentation process is outlined in Figure 4.7. From there we define boundary labels and produce volumetric tetrahedral meshes of varying resolutions. The domain boundaries are set as sketched in Figure 2.1: The basal cut corresponds to $\partial\Omega_{D}$, the epicardium to $\partial\Omega_{R}$, where the Robin boundary conditions (2.10c) are defined with a spatially varying stiffness

[TABLE]

and the endocardium to $\partial\Omega_{N}$, where we set $p_{N}(t)=p_{0}\sin^{2}(\pi t)$, representing the variation of endocardial pressure. The constants $y_{a},y_{b}$ are the vertical components of the apical and basal locations, and $\eta_{a}<\eta_{b}$ denotes the stiffness sought at the apex and base, respectively (assuming that the contact of the muscle with the aortic root is more resistant to traction than the more flexible pericardial sac and surrounding organs). In addition, since fibre and sheetlet fields for mono-ventricular geometries are not usually extracted from MRI data, we generate them using a mixed-form adaptation to the Laplace-Dirichlet rule-based method proposed in wong13 ; rossi14 .

After the S2 stimulus excites a group of cells in the lower left octant at $t=335$ ms, a spiral wave forms and sweeps around both sides of the ventricle, the two sides merging at approximately $t=415$ ms. Simultaneously, we see contraction of the apical region in the upwards direction, complemented by torsion and thickening of the ventricle wall. Figure 4.8 shows the propagation of the action potential on the deforming ventricle, with the original ventricle geometry shown with reduced opacity for comparison. The S2 stimulus occurs on the apex and the nascent scroll wave is not visible until the two arms of the wave interact. For these tests the meshsize was approximately $h=0.24\,$mm and the timestep $\Delta t=0.1\,$ms. We have employed the lowest-order finite element method $l=0$ and the stabilisation parameter is taken as $\zeta_{\text{stab}}=25$.

In addition, and similarly to the 2D case, incorporating SAD impacted the propagation of the spiral wave anisotropically. In the fibre direction, SAD led to earlier advancement of the spiral. In the transverse direction, the non-SAD case advanced earlier. Figure 4.10 shows the difference in voltage for the two cases (along with the actual voltage profile, for reference). The effect seen in the fibre direction (indicated by the white arrows) was not seen in the other directions. For these tests we have used the active strain formulation, we have included viscoelastic effects, as well as inertial contributions.

4.6 Effects due to viscoelasticity

In order to quantify the discrepancies between hyperelastic and viscoelastic effects, we conduct a series of simulations using the coupled model on a 3D slab of dimensions $50\times 50\times 10$ mm3 using a tetrahedral mesh of $h=0.25$ mm, also setting $\Delta t=0.1\,$ms, $\zeta_{\text{stab}}=25$, and ${\boldsymbol{f}_{0}}=(1,0,0)^{\tt t}$, ${\boldsymbol{s}_{0}}=(0,1,0)^{\tt t}$. These tests are conducted using the active stress formulation, and we consider inertial effects. We apply a S1 stimulus on the face $x=0$ and after $t=92\,\text{ms}$ the propagation front has reached the state shown in Figure 4.11(a), plotted on the deformed configuration (which was computed with a full electro-viscoelastic model). The boundary conditions for the viscoelasticity are of Robin type everywhere on the boundary. At that time, in panels (b,c) we depict snapshots of the approximate solutions obtained using the hyperelastic and viscoelastic models with their base-line parameter values as reported in Table 4.2, and shown over a line segment crossing the tissue slab parallel to the $x-$axis. We show profiles of the mechanical entities ($x-$components of displacement and pressure), as well as potential and $r_{3}$. For reference, we also include the results obtained using a model without SAD contributions (that is with $D_{2}=0$). We note that the curves produced without SAD are substantially lagged (as expected from the choice of diffusion parameters) with respect to the two other cases, that display no major discrepancies. The remaining panels in the figure show point-wise transients of the main mechanical and electrical fields measured on the point $(x_{0},y_{0},z_{0})=(25,25,10)$. The evolution of the electric and activation fields remains very similar in all three cases; for instance the shape of the action potential is almost not modified after adding SAD or viscous contribution and for the other fields also very subtle differences are observed (the calcium concentration was slightly shifted to the left in the hyperelastic and viscoelastic cases). The changes are more pronounced in the Frobenius norm of the Kirchhoff stress, the displacement magnitude and the pressure (panels g,h,i). These computations suggest that viscous effects will result in a decreased displacement, stress, and pressure (similar conclusions were drawn in pandolfi17 , but not in the context of models for ventricular viscoelasticity). These discrepancies, however, are qualitatively small, and this observation was robust to every parameter combination that we tested, consistent spatially and in time. The application of a viscous model also had consequences related to performance. For instance, in the tests mentioned above, the average number of Newton iterations needed to reach convergence was systematically lower in the viscous case than in the hyperelastic case. This behaviour is expected as for simple viscoelastic models the tangent problem is essentially a rescaled version of the elastic stiffness, which contributes to improving the stability of the tangent problem.

We next proceed to investigate the effects of changing the viscosity parameters. The parameter $\beta$ from (2.6) exerted minimal influence over the observed dynamics. Even for the five orders of magnitude tested, from $\beta=0.1$ ms to $\beta=10000$ ms, the differences in displacement, voltage, and all other variables were of less than 0.1%. This could be because of the low rates of change of deformation that we see in our simulations. We also tested values of $\delta$ across three orders of magnitude, from $\delta=2.26$ to $\delta=2260$ (in units N$/$cm${}^{2}\cdot$ms). As expected, increasing this quantity, thereby increasing the viscoelastic contribution to the Cauchy stress, magnified the differences between the hyperelastic and viscoelastic cases (essentially magnifying the effects seen in Figure 4.11). Additional simulations (not reported here) also showed that higher values of $\delta$ not only reduced the magnitude of $\boldsymbol{\Pi}$, $\boldsymbol{u}$, $p$, but also smoothed their profiles, reducing the distances between peaks and troughs. Even if no substantial differences were encountered in terms of conduction velocity, the calcium transients displayed generally higher values in the viscoelastic case.

4.7 Viscoelastic vs. hyperelastic effects under passive inflation and active contraction

Much more evident differences can be observed in terms of the true stress $\sigma_{f}=\mathbf{F}{\boldsymbol{f}_{0}}\cdot(\boldsymbol{\sigma}\mathbf{F}{\boldsymbol{f}_{0}})$ when plotted against the local stretch in the fibre direction, $\lambda_{f}=\sqrt{I_{4,f}}$. Such a comparison has been conducted in gultekin16 for idealised geometries, and it was specifically designed to study hysteresis effects due to viscous contributions to orthotropic passive stress. For the inflation tests we will restrict to $\beta=1$ ms and $\delta=22.6$ N$/$cm${}^{2}\cdot$ms. These values, considered in karlsen17 (and using units of [s] and [Pa s], respectively), ensure that the viscoelastic component is large enough to have a visible effect, but does not completely overwhelm the dynamics of the tissue. Here we consider the left ventricular domain used in Section 4.5 and proceed to analyse a stress-stretch response on two points near the basal surface on the endocardium and epicardium, and portrayed in Figure 4.12(a). The mechanical parameters were taken differently from those in Table 4.2; here we focus on the patient-specific constants estimated from healthy myocardial tissue at 8 mmHg end-diastolic pressure using chamber pressure-volume and strain data taken in vivo gao . The modified values for this particular test are $a=0.02096$ N$/$cm2, $b=3.243$, $a_{f}=0.30634$ N$/$cm2, $b_{f}=3.4595$, $a_{s}=0.07334$ N$/$cm2, $b_{s}=1.5473$, $a_{fs}=0.03646$ N$/$cm2, $b_{fs}=3.39$ and we set $\zeta_{\text{stab}}=10$. In the simulation we impose a sinusoidal endocardial pressure of maximal amplitude 0.1 N$/$cm2 (approximately 8 mmHg) and run a set of transient simulations over the interval from 0 to 300 ms. This configuration constitutes an inflation and deflation process where the majority of the fibres are acting in traction, whereas sheetlets work under a compression regime. Plots (b,c) in Figure 4.12 illustrate the stress-stretch response (in terms of the true stress). The behaviour on the epicardial point shows an exponential stiffening and is quite similar to what was observed in gultekin16 , as for both stress measures in the viscoelastic case there is evidence of hysteresis effects (that are, by definition, not present in the hyperelastic case). Slight deviations from the reference results in gultekin16 maybe related to the fact that we are using a full electromechanical model, a different viscoelastic contribution, and different material parameters. On the endocardial point we observe even more marked differences between the two cases, probably since we do not expect symmetry in the motion patterns for a non-ellipsoidal geometry. Other qualitative differences in the motion patterns include a more marked wall thickening, and an overall lower pressure (also more evenly spread throughout the endocardium, showing a smoother profile than the one produced in the hyperelastic case). Pressure on the epicardium was higher in the viscoelastic case.

As a final test we analyse the differences between the viscoelastic and hyperelastic case under active contraction. We employ the ventricular geometry again, imposing Robin boundary conditions on the epicardial surface and zero normal displacements on the basal cut, and use the active strain approach. We use $\zeta_{\text{stab}}=10$, and more pronounced viscoelastic effects encoded in $\beta=5$ ms and $\delta=80.25$ N$/$cm${}^{2}\cdot$ms. The active contraction of the ventricle is initiated by an ectopic beat and impose a sinusoidal endocardial pressure but with larger amplitude 0.25 N$/$cm2 and simulate the process for approximately two cycles (from 0 to 700 ms). Figure 4.13 reveals a higher asynchrony in the tissue deformation and stresses between the hyperelastic and viscoelastic case. This is observed visibly from the motion of the ventricle in panel (a), but also from the transients extracted on an epicardial point near the apex and measuring true stress in the fibre direction, the local fibre stretch, the displacement magnitude. This effect was milder in the Frobenius norm of the Kirchhoff stress tensor.

5 Concluding remarks

We have introduced a model for the active contraction of cardiac tissue. We focused on incorporating the mechanoelectric feedback through stress-assisted diffusion, accounting for a porous-media-type nonlinear diffusivity, and including inertial terms in the equations of motion. The three-field equations of motion of a viscoelastic orthotropic material are coupled with a four-variable minimal model for human ventricular action potential using both active strain and active stress approaches. We have also proposed a new stabilised mixed-primal numerical scheme written, in particular, in terms of the Kirchhoff stress. The nontrivial effects of both viscoelasticity and stress-assisted diffusion in our model suggest that they may play an important role in governing cardiac function and its response to external stimuli.

An important remark is that the active strain approach seem to be much more robust than the active stress formulation, at least in the present context. In order to obtain comparable results with the active stress approach, we had to spend quite a lot of effort finding the correct scaling in (2.3). Once this is achieved, we observed that there is no substantial difference in the output quantities. This is why our numerical tests have focused a little bit more on the active stress, since it is somewhat more challenging than the other case.

Further additions will be mostly focused on multiscale microstructural coupling, which will provide a more physiological justification of the model in terms of complex phenomena involved in mechano-electrical interactions. One example would be to include poroelastic effects representing perfusion of the myocardial tissue. Developing a thermodynamically consistent description of stress-assisted diffusion is also a pending task, in which electromechanical coupling with the surrounding torso and organs would represent another level of interaction. Such formulation under electromechanical coupling (and including nonlinear and stress-assisted diffusion), will require state-of-the-art tools of multiscale homogenisation cyron16 as well as dedicated multiscale numerical methods gandhi16 .

More confident now in obtaining accurate and reliable numerical solutions, our forthcoming contributions will target an exhaustive computational analysis of restitution curves and realistic activation patterns, e.g. accounting for Purkinje fibres and cellular heterogeneity, with the purpose of characterising spatiotemporal alternans patterns in the presence of multiple mechano-electric feedback effects. Practical applications of the present study rely on the antitachycardia pacing protocols, as well as on the (still today not completely understood) effects of mechanical loads, including cardiac massage, tissue damage, and remodelling at different scales during atrial flutter mase08 . Estimates of energy dissipation and heat production would be further investigated, widening the validity of these models to non-equilibrium thermodynamical systems. The simulation of mechanically-induced ectopic activity, as well as prediction of the risk of sudden cardiac death are also part of our long-term goals. For sake of model validation we are also interested in marrying our results to experimental observations obtained through elastography (using for instance the novel approach advanced in capilnasiu19 ).

Acknowledgements.

We express our sincere thanks to Pablo Lamata (KCL) for providing the base-line surface LV meshes, as well as Figure 4.7. In addition, fruitful discussions with Xiyao Li, Renee Miller, and David Nordsletten are gratefully acknowledged. We also thank the efforts of two anonymous referees whose many constructive comments lead to a number of improvements to the manuscript. This work has been partially supported by the Oxford Master in Mathematical Modelling and Scientific Computing, by the S19 Visiting Professorship Programme at University Campus Bio-Medico, and by the EPSRC through the research grant EP/R00207X/1.

Conflict of interest The authors declare that they have no conflict of interest.

Appendix A Details on the minimal ionic model

The ionic currents consist of three general terms, phenomenologically constructed (without particularisation to the ionic species that generate them)

[TABLE]

where the adimensional fast inward, slow inward and slow outward currents are respectively given by

[TABLE]

and the kinetics of the gating variables $\vec{r}$ are given by

[TABLE]

Here $\mathcal{H}$ is the Heaviside step function, and the time constants and steady-state values are defined as:

[TABLE]

Appendix B Mixed-primal fully discrete finite element scheme

We restrict the presentation to the active stress formulation using a smoothed model for active tension (stressing that the case of active strain and other boundary conditions follows a similar treatment). Let us denote by $\mathcal{T}_{h}$ a regular partition of $\Omega$ into simplicial elements $K$ (pair-wise disjoint triangles in 2D or tetrahedra in 3D) of maximum diameter $h_{K}$, and define the mesh size as $h:=\max\{h_{K}:K\in\mathcal{T}_{h}\}$. Let us also denote by $\mathcal{E}_{h}$ the set of interior facets of the mesh, and by $\left[\![\cdot\right]\!]_{e}$ the jump of a quantity across a given facet $e\in\mathcal{E}_{h}$. The specific finite element method we use here is based on solving the discrete weak form of the hyperelasticity equations using piecewise constant approximations of the symmetric Kirchhoff stress tensor, piecewise linear approximation of displacements, and piecewise constant approximation of (solid) pressure. The transmembrane potential in the electrophysiology equations is discretised with Lagrange finite elements (piecewise linear and continuous functions), and the remaining ionic quantities are approximated by piecewise constant functions. More precisely, we use the finite dimensional spaces $\mathbb{H}_{h}\!\subset\!\mathbb{L}^{2}_{\text{sym}}(\Omega)$, $\mathbf{V}_{h}\subset\mathbf{H}^{1}(\Omega)$, $W_{h}\subset\mathrm{H}^{1}(\Omega)$, $Q_{h}\subset\mathrm{L}^{2}(\Omega)$, $Z_{h}\subset\mathrm{L}^{2}(\Omega)^{3}$ defined (for the case of a generic-order approximation $l\geq 0$) as follows:

[TABLE]

where $\mathbb{P}_{l}(K)$ denotes the space of polynomial functions of degree $s\leq l$ defined locally on the element $K$. Assuming zero body loads, and applying a backward differentiation formula (BDF) for the time integration, we end up with the following fully-discrete nonlinear electromechanical problem, starting from the discrete initial data $v^{0}_{h},n^{0}_{h},T_{a,h}^{0}$. For each $n=0,1,\ldots,n_{\max}$: find

$(\boldsymbol{\Pi}_{h}^{n+1},\boldsymbol{u}_{h}^{n+1},p_{h}^{n+1})$ and $(v_{h}^{n+1},\vec{r}_{h}^{n+1},T_{a,h}^{n+1})$ such that

where $\zeta_{\text{stab}}$ is a positive pressure stabilisation parameter required to enforce solvability of the discrete problem. This is the tetrahedral counterpart of the finite element method for quadrilateral meshes studied in chavan07 and recently used in the context of cardiac electromechanics in ruiz18 . Notice that the boundary condition (2.10a) is incorporated as an essential condition on the displacement space, whereas the traction boundary condition (2.10b) on the remainder of the boundary $\partial\Omega_{N}$ appears naturally as the last term in the second equation of (LABEL:eq:FE).

Appendix C Linearisation of the mechanical problem

The coupling between activated mechanics and the electrophysiology solvers will be performed using a segregated fixed-point scheme. At each time step, the nonlinear algebraic sub-system for the mechanics defined by the first three equations in (LABEL:eq:FE) is linearised, adopting the following form (where the time indices have been dropped for the sake of notation).

Starting from the initial guess $(\boldsymbol{\Pi}_{h}^{k=0},\boldsymbol{u}_{h}^{k=0},p_{h}^{k=0})=(\boldsymbol{\Pi}_{h}^{n},\boldsymbol{u}_{h}^{n},p_{h}^{n})$, for $k=0,1,\ldots$ find stress, displacement and pressure increments $(\delta\boldsymbol{\Pi}_{h}^{k+1},\delta\boldsymbol{u}_{h}^{k+1},\delta p_{h}^{k+1})$ such that