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

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

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 (Augustin et al. 2016). 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 (Göktepe et al. 2013). 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 (Jaffe and Morin 2014). 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 spatiotemporal resolution to provide the required level of information. However, relatively recent studies (see, for example, Christoph et al. 2018) have used 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 multiscale and multiphysics 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, Augustin et al. 2016; Cherubini et al. 2012; Franzone et al. 2016; Costabal et al. 2017; Gizzi et al. 2015; Göktepe et al. 2013; Nobile et al. 2012; Quarteroni et al. 2017; Sundnes et al. 2014 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 (Quarteroni et al. 2017). The study of single cell and cell–cell (Lenarda et al. 2018) chemomechanical and electromechanical interactions has attempted to unveil some of the underlying complex features of the cardiac function, and different multifield nonlinear models have been gradually generalising classical approaches such as the monodomain equations and Fick’s law of diffusion. In particular, fractional diffusion (Cusimano et al. 2018), nonlinear diffusion (Hurtado et al. 2016), and stress-assisted diffusion formulations (Cherubini et al. 2017) 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 (Phadumdeo and Weinberg 2018). At the same time, macroscale incompressibility, orthotropic, and hysteretic mechanical features have been shown to fully characterise the human cardiac tissue under multiaxial loading tests (see, for example, Gültekin et al. 2016 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 characteristics describing the viscous response due to matrix, fibre, sheet, and fibre–sheet couplings through four dedicated dashpots (Gültekin et al. 2016). Also in this case, a porous medium motivation has been advanced in Yao et al. (2012), 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 (Holzapfel and Gasser 2001), motivated by a rheological analogue from Simo 1987, and endowed with equilibrium and non-equilibrium contributions (Lubliner 1985) 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 (Klepach and Zohdi 2014), but has only recently been adapted for active biological media undergoing reaction–diffusion excitation (Cherubini et al. 2017), and more tailored for cardiac models in Loppini et al. (2018). 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 Ruiz-Baier (2015) and used more recently in Cherubini et al. (2017), Garcia-Blanco et al. (2019), and Ruiz-Baier et al. (2019). In the present case, our model and our numerical method include a three-field elasticity formulation (variationally based on a modification of the Hu–Washizu principle (Lamichhane et al. 2006)) that states the governing equations in terms of stress–displacement–pressure, motivated by the desire to avoid 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 Chavan et al. (2007) 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 (Bueno-Orovio et al. 2008). 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 (Rossi et al. 2014), the active stress approach is somewhat simpler and more naturally incorporated in already existing models for passive deformation (Giantesio et al. 2019). 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 Sect. 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 Sect. 4. We close with a summary and some remarks on model limitations and ongoing extensions, collected in Sect. 5.

Mathematical model

Finite-strain cardiac mechanics

Let , , denote a deformable body with a piecewise smooth boundary , considered in its reference configuration, and let denote the outward unit normal vector on . The kinematical description of finite deformations regarded on a time interval is made precise as follows. A material point in is denoted by , whereas will denote the displacement field defining its new position in the deformed configuration. The tensor is the gradient (applied with respect to the fixed material coordinates) of the deformation map; its determinant, denoted by , measures the solid volume change during the deformation; and and are, respectively, the right and left Cauchy–Green deformation tensors on which all strain measures will be based (here the superscript denotes the transpose operator). The first isotropic invariant ruling deviatoric effects is the scalar , and for generic unitary vectors , the scalars , , are pseudo-invariants of measuring direction-specific stretch (Ciarlet 1988).

The triplet represents an orthogonal coordinate system pointing in the local direction of the muscular cardiac fibres, transversal sheetlet compound, and normal cross-fibre direction . Note that the system is restricted to 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 Holzapfel and Ogden (2009), 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, respectivelywhere a, b are material constants associated with the isotropic matrix response, and rule the directional behaviour of the material along myocardial fibres, and account for the cross-contribution of the fibre–sheet directions, and encapsulate the shear effects in the fibre–sheet plane. Moreover, the field p denotes the solid hydrostatic pressure, and we use the notation if or zero otherwise, for a generic real-valued function u. This modelling choice is appropriate given 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 (Pezzuto et al. 2014). We remark here that the particular mechanisms of soft tissue anisotropic mechanical behaviour are still under investigation (Humphrey et al. 2014), and the chosen formulation may not be the most general one. 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 (McEvoy et al. 2018).

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, Sundnes et al. 2014): we assume that the first Piola-Kirchhoff stress tensor decomposes aswhere the active stress component acts differently on each local direction with an intensity depending on the scalar field of active tension , that synthesises (in an homogenised sense) the biochemical state of myocytes (and whose dynamic behaviour will be specified later on). Then,where are positive constants representing the variation of activation on each specific direction, as proposed in Dorri et al. (2006), and are the fibre, sheetlet, and cross-fibre stretches. Setting appropriate models for 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 (Pezzuto et al. 2014). Details of other anisotropic activation forms can be found, for instance, in Rossi et al. (2014) for active strain and in (Usyk et al. 2000, 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 since a stress component on this direction would counteract wall thickening mechanisms (Dorri et al. 2006). Moreover, the intensity of the active tension effect on the cross-fibre direction is assumed to be substantially smaller than that appearing on the off-diagonal component ; see Table 2. Also note that some references do not include a normalisation with local stretches in each term of . Finally, these constitutive laws are usually not derived from a thermodynamic potential.

Table 2

Model parameters for the electro-viscoelastic model (2.8), (2.11), (2.16), (2.14). Values are taken from Cherubini et al. (2017), Gao et al. (2015), Rossi et al. (2014), Bueno-Orovio et al. (2008), and the transmembrane potential v is dimensionless

Viscoelasticity constants
a=0.236	[N/cm^2]	a_f=1.160	[N/cm^2]	a_s=3.724	[N/cm^2]	a_fs=4.010	[N/cm^2]
b=10.81	[–]	b_f=14.15	[–]	b_s=5.165	[–]	b_fs=11.60	[–]
p_0=0.1	[N/cm^2]	β=10	[ms]	δ=22.6	[N/cm^2 ms]	ζ_stab=0.25	[–]
η_a=0.001	[N/cm^2]	η_b=0.01	[N/cm^2]	κ_sn=0.6	[–]	κ_nn=0.03	[–]
ρ_0=0.001	[N/cm^2]

Active strain model. Next we recall the active strain model for ventricular electromechanics (see, for example, Cherubini et al. 2008). 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, , withThe coefficients , with , 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 potentially 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 , which determines a new (and fictitious, or virtual) intermediate configuration that is regarded as a reference for the elastic deformation (Pezzuto et al. 2014). Therefore, the strain energy function and the first Piola-Kirchhoff stress tensor (after applying the active strain decomposition) are functions of only, and read, respectivelyAs 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 (Pezzuto et al. 2014) (these will be the terms appearing on the second diagonal block of the weak formulation from Sect. 3, the block corresponding to displacements); however, the overall problem will remain of a saddle-point structure. The modified elastic invariants are functions of the coefficients , as well as of the invariant and pseudo-invariants in the following manner (Rossi et al. 2012)Such dependencies are a consequence of assuming isochoric active deformations (Pezzuto et al. 2014), i.e. , justified by the fact that the volume of the cardiomyocytes does not vary substantially during contraction. Besides, following Rossi et al. 2012, previous expressions are obtained by assuming and making use of the fact that , as well as that , with Accordingly, the active strain, and consequently the force associated with 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 (Rossi et al. 2014).

Viscoelasticity and equations of motion

Extension and shear tests demonstrate the importance of incorporating viscoelastic effects in models for cardiac passive mechanics (Gültekin et al. 2016). In the heart, the extracellular fluid filtrating through the elastic solid is one of the main generators of the viscoelastic effects of the tissue (Yao et al. 2012). Viscous effects are also tied to cross-bridge processes identified in Ca activated fibres (Maughan et al. 1998) 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 (Holzapfel and Gasser 2001). From the viewpoint of kinematics, it suffices to relate stress to strain rates. Decomposition of the spatial velocity gradient into the rate of deformation and spin tensors yields the relationand a simplified rheological Kelvin–Voigt model for the viscous component of the Cauchy stress can be defined as follows (see, for example, Karlsen 2017):which depends on the history of the isotropic contribution to the Cauchy stress. Here, are model parameters. In this way, after a pull-back operation, we see thatis the total first Piola-Kirchhoff stress tensor that includes defined from either (2.1)–(2.2) or (2.5), and the viscoelastic contributions.

More advanced rheologies can be easily incorporated in the context of active stress formulations as done in, e.g. (Katsnelson et al. 2004), as the generalised Hill–Maxwell model recently proposed in Cansiz et al. 2017, as in the perturbed equations of harmonic wave motion using springpot-based models with fractional order derivatives from Capilnasiu et al. 2019, or as in the thermodynamical electro-viscoelastic models that use statistical fibre distributions (Pandolfi et al. 2017; Gizzi et al. 2018). 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, where are the reference and current medium density, is a smooth vector field of imposed body loads, denotes the second time derivative, and the divergence operator in (2.8a) acts on the tensor fields row by row. The balance of angular momentum translates into the condition that the Kirchhoff stress tensor must be symmetric, which is in turn encapsulated into the momentum and constitutive relations (2.8a), (2.1), (2.5).

Definingas the contribution to the Kirchhoff stress that does not involve pressure, we then haveStating 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, Eqs. (2.8a)–(2.8b)–(2.9) will be supplemented with mixed normal displacement and traction boundary conditions where , , 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 . The term in (2.10b) denotes a (possibly time dependent) prescribed boundary pressure associated with endocardial blood pressure, which is uniform over the deformed counterpart of , and it is applied in the normal direction to the epicardium in the deformed configuration. However, owing to Nanson’s formula (Ciarlet 1988), 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 . 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 ). More sophisticated boundary conditions that consider an interaction with the pericardium can be also imposed (Fritz et al. 2014). A sketch of a mono-ventricular domain specifying boundary surfaces and fibre directions is depicted in Fig. 1.

Fig. 1

Schematic representation of a mono-ventricular domain where (2.10a) is imposed on the basal cut, (2.10b) on the endocardial surface, and (2.10c) on the epicardium. The left panel depicts the fibre field and the right panel the sheetlet directions (in this case, parallel to the normal direction of the epicardium)

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, for example, Franzone et al. 2014), 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 is encoded in the vector (here we use instead of bold to denote vector fields of dimension other than d) Here, is the ratio of membrane area per tissue volume, and is a spatiotemporal external stimulus applied to the medium. We will adopt the minimal model for human ventricular action potential, proposed in Bueno-Orovio et al. (2008) 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 Gizzi et al. (2013). That model was, however, tailored originally for the case of isotropic conductivity , 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 (Cherubini et al. 2017; Loppini et al. 2018).

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

Boundary and initial conditions for (2.11) correspond to 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 Bueno-Orovio et al. (2008), 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 2), assuming that it is consistent throughout the cardiac wall. Extension to the heterogeneous case can be readily incorporated into the present framework.

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 Cherubini et al. (2017) and later used for simplified 2D cardiac electromechanics in Loppini et al. (2018)). 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 constitutes a strain-based modification of tissue conductivity, also referred to as geometric feedback in Franzone et al. (2016))where the nonlinear conductivity (self-diffusion depending on v) accounts for porous media electrophysiology following the development in Hurtado et al. (2016), but appropriately modified to incorporate information about preferred directions of diffusivity according to the microstructure of the tissue (encoded in the second term defining ). The parameter signifies the usual diffusion for isotropic materials, whereas and 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. , as in Gizzi et al. (2017), Ruiz-Baier et al. (2019)) eventually leads to very slight modifications in conduction patterns, and we have therefore decided not to include it. Tuning 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 (Aifantis 1980) 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 (Lenarda et al. 2018). 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 (Ly and Weinberg 2018), as well as micro-invaginations on the cell membrane known as microtubules and microdomains (Miragoli et al. 2016). 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, as well as validated using an experimental dataset.

It is important to remark that the solvability of the monodomain equations (2.11a)–(2.11b) depends on the properties of . 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 Cherubini et al. (2017). Here, we perform a much lighter calibration, as mentioned later in Sect. 4. Comparisons between the effects of SAD and the more conventional mechanoelectrical feedback through stretch-activated currents have been reported in Loppini et al. (2018).

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 formwith , and , where are positive model constants. As calcium concentration is not readily available in the phenomenological cellular model we are employing, we use as a proxy for intracellular calcium (Bueno-Orovio et al. 2008). In addition, a linear dependence on the calcium proxy and on the local stretch is 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 in terms of the microscopic cell shortening is expressed as follows (see, for example, Barbarotta et al. 2018):and the specific relation between the myocyte shortening and the dynamics of slow ionic quantities (in the context of our phenomenological model, only ) is made precise using the lawwhich does not require an explicit dependence on local fibre stretch, as the sliding of myofilaments is driving the dynamics of the functions . We employ the nonlinear reaction term , and we make the distinction that and characterise the evolution of the activation in the approaches of active stress and active strain, respectively.

Numerical method and implementation

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 ) 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 , find as well as such thatwhere , and where the case for an active stress formulation necessitating an active tension model is addressed similarly (however, the regularity of is then ). Theoretical aspects regarding the coupling of elasticity and stress-assisted diffusion problems have been recently addressed in the context of mixed-primal and mixed-mixed formulations in Gatica et al. (2018), 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 (3.1). Details on the finite-dimensional spaces and linearisation procedure are laid out in Appendix 2.

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 (Lamichhane and Stephan 2012); 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) (Chavan et al. 2007). 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 , the matrix system associated with (B.1) 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 (Loppini et al. 2018)]. The importance of casting the equations of motion in terms of the coupling variables has been already emphasised in Ruiz-Baier (2015) 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.

Solver structure and implementation details

According to the separation (through an outer fixed-point scheme) 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 and 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 Sect. 4.1, which will be solved with the direct method SuperLU), and the iterates are terminated once a tolerance of (imposed on the 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 (Quarteroni et al. 2017). All routines have been implemented using the finite element library FEniCS (Alnæs et al. 2015).

Computational results

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 Sect. 3. Let us consider the following closed-form solutions to a steady-state counterpart of the variational form (3.1) for the electromechanics equations, also assuming the absence of viscoelastic effects, and defined on the domain with the fibres/sheetlets defined as Then, the Kirchhoff stress and 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 equationsNote also that the incompressibility constraint for this test is , where 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 and . Kirchhoff stress and pressure errors are measured in the norm, whereas for the remaining variables, the errors are measured in the norm. The obtained error history is reported in Table 1, where we observe an asymptotic 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 (Chavan et al. 2007).

Table 1

Test 1: Error history (errors on a sequence of successively refined grids and convergence rates) associated with the mixed finite element method (B.1) applied to a steady-state electromechanical coupling under active stress, and using different polynomial degrees

(a) Hyperelasticity variables
DoF	h	‖Π-Π_h‖_0,Ω	rate	‖u-u_h‖_1,Ω	rate	‖p-p_h‖_0,Ω	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

Parameter calibration

For the following 2D simulations, we will initially consider tissue slabs of mm, and set fibre and sheetlet directions simply as , . The initiation, maintenance, prevention, and treatment of so-called reentrant waves are major focus of the current research due to their implication in atrial and ventricular fibrillations (Franzone et al. 2014). 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 non-dimensional amplitude 3 and duration 3 ms. The S1 stimulus is a planar wave created by exciting the entire left edge of the tissue, and the S2 stimulus is a square wave created by exciting the bottom left quadrant at t = 330, ms. In 3D, the same general protocol can be used, but the S1 stimulus excites the entire bottom section (below some value of the coordinate) while the S2 stimulus excites the bottom left octant at ms. 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 Fig. 2. 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.

Fig. 2

Evolution of voltage after S2 stimulus, of S1-S2 protocol, showing formation of a reentrant spiral wave on the deforming viscoelastic tissue, computed using the active stress approach

Next, since we are using the active stress formulation in this case, we proceed to evaluate , the parameters governing active tension in (2.14), and , the stiffness parameter from (2.10c). We conduct a simple sensitivity analysis by increasing or decreasing either , or by one order of magnitude, holding the others constant at their reference values (, , and N/cm, as listed in Table 2). This simple analysis therefore does not test for compounding or interaction effects. We also consider a smaller slab of size mm. The parameter contributes to producing smoother active tension profiles, while controls the range of their magnitude. These effects are visible in Fig. 3. We found that larger values of produced smoother gradients in pressure and stress, while larger values of produced, in average, higher magnitude displacement, Kirchhoff stress, and pressure, as well as some more subtle changes in ionic quantities. Parameter determines the stiffness of the springs supporting the tissue, and so decreasing 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 (E−4 to ). The effects on ionic entities were even smaller, for both the hyperelastic and viscoelastic cases, and therefore, plots are not shown.

Fig. 3

Profiles of taken across a smaller slab of tissue at and . These plots evaluate the effect of and

Computational experiments reveal a window of values of for which our method converges. In the 2D hyperelastic case, we found that the upper bound for is approximately E−2 cm(s Pa), with the linear solver failing to converge for larger values. In these simulations, the Kirchhoff stress achieved an norm of between 0.006 and 0.6. In turn, the viscoelastic case was able to accept slightly larger values of , up to E−2 cm(s Pa), with the 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 Cherubini et al. (2017).

The numerical method used for these tests is characterised by the time step, mesh size, polynomial degree, and stabilisation constant  ms,  mm, , , respectively.

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 (Chavan et al. 2007, Test II)), where we set an active tension of . This test involves applying an upward in-plane shear load 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 (see the sketch in Fig. 4a), and the fibre and sheetlet fields are and , respectively. Secondly, we consider a 3D system suggested in (Land et al. 2015, Test I) as a simple benchmark for passive cardiac mechanics, and therefore, we set . The problem consists in computing the deformation of a point at the right end of a beam defined by the domain   mm (see Fig. 4b), where the fibre direction is . Instead of (2.1), the material is characterised by the transversally isotropic strain energy density proposed by Guccione et al. (1995) [which is the material law used in the benchmark test from Land et al. (2015)]: , with , where  kPa, , , , and the denote entries of the Green–Lagrange strain tensor , rotated with respect to a local coordinate system aligned with . The beam is clamped at the face , a pressure of  kPa is imposed on the bottom face , 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.

Fig. 4

Domain sketches and sample meshes for the deflection of Cook’s membrane for an Holzapfel–Ogden material with constant active stress (a) and deflection of a 3D beam for a Guccione–Costa–McCulloch material with the active stress component set to zero (b)

The outcome of these tests in Fig. 5a, 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 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 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 Land et al. (2015)). Moreover, Fig. 5c, d shows 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, , and the observed behaviour also confirms that our method is locking-free.

Fig. 5

Convergence of the deflection of Cook’s membrane for an Holzapfel–Ogden material with constant active stress (a, c) and deflection of a 3D beam for a passive Guccione–Costa–McCulloch material (b, d). Maximal vertical deflection with respect to the mesh resolution for different numerical schemes (a, b), and different values of the stabilisation constant (c, d)

Stress-assisted diffusion and conduction velocity

In addition to determining a suitable parameter range for that ensures solvability of the discrete monodomain equations, we also investigated the effect of on the tissue’s response to spiral wave dynamics. As in the second part of Sect. 4.2, this time the domain is a square slab of width 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 and the sheetlets , and the Holzapfel–Ogden material law is considered. The mesh size is approximately  mm and the time step is ms. We use the lowest-order finite element method and the stabilisation parameter is .

Figure 6 shows the differences in the ionic quantities between simulations with a very small contribution of SAD ( 1 E−5 m(s Pa)) and a more prominent, but still mild SAD contribution (7.5E−3 cm(s Pa)). The snapshots correspond to the time 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 Fig.  7a, b, we see that conduction velocity was higher for larger values of (meaning a larger SAD contribution). When the wave first emerged, the peak action potential was more advanced for the case of reduced , but the large 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 Fig.  7a, b revealed that the overall profiles of the ionic quantities were highly similar between the two cases compared in Fig. 6, but differed in the speed at which they are transported through the tissue.

Fig. 6

Differences in ionic quantities from varying SAD parameter at . Quantities indicate the profiles with 7.5E−3, and the profiles associated with 1.0E−5 (dimensions are as in Table 2)

Fig. 7

a, b Propagation of action potential v and active tension , measured by taking the profile over a horizontal line segment crossing the upper half of the tissue at . Comparison is provided for two different values of . c, d, e Effect of on the potential wave at in the viscoelastic case

We also remark that the effect of changing conduction velocities was not spatially consistent. For some parameter values, 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. Figure 7c–e shows a comparison of the spiral wave in the viscoelastic case for three different values of . The upper right area of the spiral is slightly more vertical in the simulation with a larger value of than in the other cases, suggesting that propagation of the voltage is somewhat hindered in that 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 spatiotemporal discretisation. In Table 3, we include the results of a simple convergence test for conduction velocity, similar to the benchmark test conducted in Ruiz-Baier et al. (2019). 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 Ruiz-Baier et al. (2019), 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 (, 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.

Table 3

Convergence of conduction velocity with respect to temporal and spatial discretisation

Convergence of conduction velocity, mm/ms
DoF	h (mm)	Δt=0.3ms	0.1ms	0.05ms	0.01ms
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

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 Warriner et al. (2018), Lamata (2018)] and rescale them using approximately the same dimensions as the idealised ventricles studied in Ruiz-Baier et al. (2019). The segmentation process is outlined in Fig. 8. From there we define boundary labels and produce volumetric tetrahedral meshes of varying resolutions. The domain boundaries are set as sketched in Fig. 1: the basal cut corresponds to , the epicardium to , where the Robin boundary conditions (2.10c) are defined with a spatially varying stiffnessand the endocardium to , where we set , representing the variation of endocardial pressure. The constants are the vertical components of the apical and basal locations, and 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 Wong and Kuhl (2014), Rossi et al. (2014).

Fig. 8

Segmentation and mesh personalisation process from Warriner et al. (2018), Lamata (2018). Semi-automatic segmentation by 3D extrapolation (yellow surface and contours) of 2D segmentation contours (red contours and projections) (a); surface mesh template (b); and resulting mesh (white surface) overlaid with the segmentation surface colour coded by the distance between them (jet colour map, from 0 mm in blue to 1 mm in red) (c). Used with permission

After the S2 stimulus excites a group of cells in the lower left octant at ms, a spiral wave forms and sweeps around both sides of the ventricle, the two sides merging at approximately ms. Simultaneously, we see contraction of the apical region in the upwards direction, complemented by torsion and thickening of the ventricle wall. Figure 9 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 mesh size was approximately  mm and the time step ms. We have employed the lowest-order finite element method and the stabilisation parameter is taken as  (Fig. 10).

Fig. 9

Evolution of voltage after S2 stimulus (at ), showing formation of a scroll wave on a contracting ventricle, using the active strain model. The shadow of the undeformed ventricle geometry is shown for comparison

Fig. 10

Snapshot at of field variables plotted on the deformed domain and less opaque undeformed mesh. Here we have also used the active strain approach

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 11 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 and we have included viscoelastic effects, as well as inertial contributions.

Fig. 11

Effect of SAD on spiral wave propagation, using the active strain formulation. Panels a, b show voltage and c shows the difference between the SAD and non-SAD cases (which has a different scale). The action potential wave using SAD moved along the fibre direction ahead of the non-SAD case

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 mm using a tetrahedral mesh of  mm, also setting  ms, , and , . These tests are conducted using the active stress formulation, and we consider inertial effects. We apply a S1 stimulus on the face , and after , the propagation front has reached the state shown in Fig. 12a, 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 baseline parameter values as reported in Table 2, and shown over a line segment crossing the tissue slab parallel to the axis. We show profiles of the mechanical entities (components of displacement and pressure), as well as potential and . For reference, we also include the results obtained using a model without SAD contributions (that is with ). 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 pointwise transients of the main mechanical and electrical fields measured on the point . 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 Pandolfi et al. (2017), 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 linearised system.

Fig. 12

Comparison of field variables between hyperelastic and viscoelastic cases on a line parallel to the x axis (sketched in a) taken at (b, c); and pointwise evolution of field variables on the point (d–i) for the cases of hyperelasticity without SAD, with the baseline case of SAD but without viscous stresses, and the viscoelastic case (line, dashed, and dashed-dotted curves, respectively). For these tests, we have used the active stress formulation

We next proceed to investigate the effects of changing the viscosity parameters. The parameter from (2.6) exerted minimal influence over the observed dynamics. Even for the five orders of magnitude tested, from  ms to  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 across three orders of magnitude, from to (in units N/cmms). 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 Fig. 12). Additional simulations (not reported here) also showed that higher values of not only reduced the magnitude of , , 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.

Viscoelastic versus hyperelastic effects under passive inflation and active contraction

Much more evident differences can be observed in terms of the true stress when plotted against the local stretch in the fibre direction, . Such a comparison has been conducted in Gültekin et al. (2016) 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  ms and  N/cmms. These values, considered in Karlsen (2017) (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 Sect. 4.5 and proceed to analyse a stress–stretch response on two points near the basal surface on the endocardium and epicardium, portrayed in Fig. 13a. The mechanical parameters were taken differently from those in Table 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 et al. 2015). The modified values for this particular test are  N/cm, ,  N/cm, ,  N/cm, ,  N/cm, , and we set . In the simulation, we impose a sinusoidal endocardial pressure of maximal amplitude 0.1 N/cm (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 Fig. 13 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 Gültekin et al. (2016), 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 Gültekin et al. (2016) may be 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.

Fig. 13

Comparison between hyperelastic and viscoelastic deformation under passive inflation. True stress in the fibre direction , measured according to local stretch on two points on the epicardium (b) and endocardium (c) [points indicated in panel (a)]. The plots in panels (d, e) show transients of mechanical outputs (Frobenius norm of the Kirchhoff stress, true stress on fibre direction, local stretch, and displacement magnitude) at the point ; and plots (f, g) display their counterparts in point . For these tests we have used only inertial effects and passive hyperelastic or viscoelastic contributions

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 , and more pronounced viscoelastic effects encoded in  ms and  N/cm2 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/cm and simulate the process for approximately two cycles (from 0 to 700 ms). Figure 14 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.

Fig. 14

Comparison between hyperelastic and viscoelastic effects under contraction, using the active strain approach and  ms and  N/cm2 ms. Panel a has side views of the deformed domain for viscoelastic (hollow blue) and hyperelastic (dark red dots) at three different times, and panels b, c show transient of mechanical outputs extracted from a point on the lower epicardial surface

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 non-trivial 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 seems 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 mechanoelectrical 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 (Cyron et al. 2016) as well as dedicated multiscale numerical methods (Gandhi and Roth 2016).

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 (Gizzi et al. 2013) in the presence of multiple mechanoelectric 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 (Masé et al. 2008). 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, is also part of our long-term goals. For the 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 Capilnasiu et al. (2019).