We will assume that the body is a hyperelastic material and its motion will be described using finite kinematics. We will adopt an indicial notation where repeated indices indicate summation. We identify the relationship between material (reference), X_I, and spatial (deformed), x_i, coordinates via the smooth map x_i(X_I). The deformation gradient tensor F_iI = ∂x_i/∂X_I allows to determine further properties of the continuum's motion. We indicate with J = detF_iI the Jacobian of the map and with C_IJ = F_kIF_kJ and B_ij = F_iKF_jK the right and left Cauchy-Green deformation tensors, respectively. We assume that the generic myocardial fiber direction (the unit vector characterizing the microstructural property of the continuum body) in the material configuration, a_I, is mapped to the deformed configuration as a_i = F_iJa_J such that we can define the current fiber a_i = a_I/λ. Following the standard frame indifference mechanical framework (Spencer, 1989), these quantities are related to the invariants of the deformation in the following manner

The principal invariants I₁ and I₂ rule the deviatoric response of the medium, the third invariant I₃ quantifies volumetric changes of the material, while the fourth pseudo-invariant I₄ measures the directional fiber stretch, λ. This last entity is intrinsically directional, so for two-dimensional models, we will simply assign a horizontal myocardial direction (1, 0)^T. In what follows, the symbol δ_ij denotes the second-order identity tensor.

As anticipated above, we will base our model on the stress-assisted diffusion formulation from Cherubini et al. (2017). We do however, generalize the governing equations adopting a more accurate nondimensional three-variable model of cardiac action potential (AP) propagation introduced in Fenton and Karma (1998b), and we will account for SAC (Panfilov and Keldermann, 2005), that were not considered in Cherubini et al. (2017). Even though several more physiological assumptions could be made, here we will focus on a purely phenomenological approach.

In the deformed configuration, the electrophysiological model consists of three variables: the membrane potential u, and a fast and slow transmembrane ionic gates v, w. They satisfy the following RD system

where Neumann zero-flux boundary conditions are imposed for Equation (1a), i.e., [d_ij∂u/∂x_j]n_i = 0, where n_i is the outward normal on the domain boundary. System (1) describes the propagation of a normalized dimensionless membrane potential, which can be mapped to physical quantities as u = (V_m−V_o)/(V_fi−V_o) (see Fenton and Karma, 1998b for details as modified Beeler-Reuter fit) where V_m stands for the physical transmembrane potential, V_o is the resting membrane potential and V_fi represents the Nernst potential of the fast inward current. In Equation (1a), the total transmembrane density current, I_ion(u, v, w), is the sum of a fast inward depolarizing current, I_fi(u, v), a slow rectifying outward current, I_so(u), and a slow inward current, I_si(u, w), given by

where τv-(u)=Hvτv1-+(1-Hv)τv2- is the time constant governing the reactivation of the fast inward current, and H_x = H_x(u−u_x) is the standard Heaviside step function. I_ext is the space and time-dependent external stimulation current with amplitude Iextmax. All model parameters are collected in Table 2.

The mechanical problem, stated also on the current configuration and occupying the domain Ω(t), respects the balance of linear momentum and mass, written in terms of displacement, φ, and pressure, p, and set in a quasi-static form. The problem is complemented with displacement and traction boundary conditions set on two different parts of the boundary Γ_D or Γ_N:

where ρ₀, ρ and dV^,dv^ are the densities and volumes of the solid in the undeformed and deformed configurations, respectively. In Equation (3b), φ~(t) is a known (possibly time-dependent) displacement and in Equation (3c), t~i(t) is a (possibly time-dependent) traction force. In both cases, the tissue is stretched up to a maximum level of 20% of the resting length such to activate all MEF components. In addition, the time-variation of the imposed boundary conditions is much slower than the governing dynamic physical processes, and therefore a quasi-static mechanical equilibrium is maintained.

The two sub-problems (Equations 2, 3) are completed via the following mixed constitutive prescriptions for incompressible isotropic hyperelastic materials (J = 1):

Equation (4a) specifies a constitutive form for the Cauchy stress tensor (total equilibrium stress in the current deformed configuration) highlighting two multiscale contributions on the tissue deformation. First, the passive material response follows that of an incompressible Mooney-Rivlin hyperelastic solid and it is characterized by two stiffness parameters c₁ and c₂; and secondly, the active component contributing to the total stress in the form of an additional hydrostatic force with amplitude T_a. The dynamics of T_a are described by Equation (4a), where the constant k_Ta modulates the amplitude of the active stress contribution, while ϵ(u) is a contraction switch function: ϵ(u) = ϵ₀ if u < 0.005, and ϵ(u) = 10ϵ₀ if u ≥ 0.005.

Equation (4c) characterizes the stress-assisted diffusion contribution describing the effect of tissue deformation on the AP spreading. The parameter D₀ represents the usual diffusion coefficient for isotropic media, i.e., diffusivity = [L² T⁻¹], while D₁ and D₂ introduce the impact of mechanical stress through linear and nonlinear contributions, respectively, on the diffusive flux. Accordingly, D₁ and D₂ have units of [L² T⁻¹ P⁻¹] and [L² T⁻¹ P⁻²], respectively. We also remark that Equation (4c) reduces to the characterization of the classical diffusion equation for D₁ ≡ D₂ = 0.

Finally, Equation (4d) describes the stretch-activated current contribution (which is usually adopted as the sole MEF effect). The term Isac(λ, u) affects the ionic (reaction) currents in the electrophysiological system and is formulated as a linear function of the membrane potential u and the fiber stretch λ. Here, G_s modulates the amplitude of the current, usac represents a referential (resting) potential while, H_sac is a switch activating this additional reaction current only when the myocardial fiber is elongated, i.e., H_sac = 1 for λ ≥ 1 and H_sac = 0 for λ < 1.

We also introduce the definition of spiral tip (core of the spiral wave) as the point with instantaneous null velocity (see Fenton and Karma, 1998b for details). In practice, for two-dimensional domains, we choose an isopotential line of constant membrane voltage, u(R_I, t) = u_iso, where R_I = x_tipX_I + y_tipY_I represents the position vector in the reference undeformed configuration identifying the boundary between depolarized and repolarized regions. Accordingly, the spiral tip can be defined as the point in space where the excitation front meets the repolarization waveback of the action potential, conforming with the operative definition:

We numerically identify the tip coordinates (x_tip, y_tip) by considering u_iso = 0.5 with tolerance of 10⁻⁴.

The electromechanical problem is rewritten in the undeformed configuration and subsequently computationally solved via a finite element method. Even if the model originates as an extension of our contribution in Cherubini et al. (2017), the numerical method employed here is simpler, as we do not solve for stresses explicitly but rather postprocess them from the computed discrete displacements. The overall numerical scheme for active stress electromechanics with SAC is therefore not precisely novel, but we will still provide a few details for sake of completeness of the presentation and future reproducibility of results. Further details could be found in e.g., Ruiz-Baier (2015). We discretize displacements with vectorial piecewise quadratic and continuous polynomials, and the pressure field using piecewise linear and discontinuous elements. All remaining unknowns (associated to the electrophysiology and to the active tension) are approximated using piecewise linear and continuous elements. Let us then consider a regular, quasi-uniform partition Th of Ω(0)¯ into triangles T of diameter h_T, where h=max{hT:T∈Th} is the meshsize. The finite element spaces mentioned above are defined as (see e.g., Quarteroni and Valli, 1994)

for the case of clamped boundaries at Γ_D(0).

Let us also construct an equispaced partition of the time domain 0=t0<t1=Δt<⋯<tM=tmax. The coupled problem is solved sequentially between the mechanical and electrochemical blocks. A description of the needed computations at each time step tⁿ is as follows:

Step 1: From the known values uhn,vhn,whn,Ta,hn,Dhn,λhn, find uhn+1,vhn+1,whn+1,Ta,hn+1 such that

for all (ψhu,ψhv,ψhw,ψhTa)∈[Vh]4. This scheme for the electric/activation system is given in a first-order semi-implicit form: the nonlinear reaction terms and the coupling stress-assisted diffusion are taken explicitly, while the linear part of diffusion is advanced implicitly. Here

are the explicit approximation of the stress-assisted diffusivity and of the stretch in the fiber direction, all in the reference configuration.

Step 2: Given the activation value Ta,hn+1 computed in Step 1 of this iteration, solve the nonlinear elasticity equations

where

is the second Piola-Kirchhoff stress tensor.

Step 3: The solution of the problem in Step 2 uses a Newton-Raphson method whose iterations are terminated once the energy residual drops below the relative tolerance of 10⁻⁶. The solution to each linear tangent problem is conducted with the BiCGSTAB method preconditioned with an incomplete LU(0) factorization. The iterations of the Krylov solver are terminated after reaching the absolute tolerance 10⁻⁵. The residual computation for the mechanical problem also contains the terms arising from time-dependent displacement or traction boundary conditions, which also need to be assigned at each timestep. For instance, in an uniaxial test (denoted dynamic displacement in the examples below), the left segment of the boundary is clamped (zero displacements are imposed), the bottom and top edges are subject to zero normal stress, and the right edge is pulled according to the displacement φ~(t)=[0.2Lsin2(π/400t),0]T.

All tests are conducted using a two-dimensional slab of dimensions L × L = 6.2 × 6.2cm², which is the same configuration used to produce the dynamics analyzed in Fenton and Karma (1998b). The computational domain is discretized with a structured triangular mesh of 10,000 elements. After a mesh convergence test involving conduction velocities and reproducing the expected values for planar excitation waves reported in Fenton and Karma (1998b), we proceeded to fix the temporal and spatial resolutions to Δt = 0.1ms, h = 0.062cm, respectively. A representative example of the mesh is provided in Figure 4, plotted in the deformed configuration under both traction and displacement boundary conditions and highlighting the spiral wave resolution. All numerical tests were carried out using the open-source finite element library FEniCS (Alnæs et al., 2015).

We start analyzing the parameter space associated to the two MEF contributions in our model. That is, the stress-assisted coefficients D₁, D₂ and the SAC amplitude G_s. The study will be restricted to a static homogeneous stretched state (e.g., a uniaxial Dirichlet boundary condition φ = [0.2L, 0]^T set on the right edge of the domain). All remaining material and electrophysiology parameters will be kept constant, except that we fix the relative influence of the nonlinear contribution in the stress-assisted diffusion, by setting D₂ to be one order of magnitude smaller than D₁. This configuration will highlight MEF effects in a minimal, but still comprehensive manner.

Figure 5 portrays the conduction velocity obtained for all combinations of (D₁, G_s) on the parameter space. The quantity is measured as the wave-front velocity of a planar excitation wave along its propagation. The plot illustrates the variability of the recorded CV amplitude (in the range 0.25–0.5 m/s) according to the MEF coupling intensity variation and to histogram measures in Figure 2. In particular, starting from a physiological baseline of 0.42 m/s, when neither SAC nor SAD is present (D₁ = 0, G_s = 0), we observe a net increase of CV for (D₁ = 0, G_s > 0) while we recover CV decrements for (D₁ < 0, G_s = 0). This specific aspect reproduces what is expected from experimental evidence, i.e., MEF decreases the CV of the excitation wave (Ravelli, 2003).

Besides, for higher values of G_s, we obtain two unexpected results. First, for G_s > 0.15 we observe a decrement of CV for different values of D₁. Second, for the particular combination (D1<-10-4,Gs>0.15) the wave disappears from the domain or annihilates due to excessive activation (see e.g., side panels in Figure 5 or the top row in Figure 8). Consequently, we are not able to measure any propagation (which reflects in the combinations with × of the figure). This last result is somehow counterintuitive since, as evidenced by Figure 1, we experimentally experience a complete depolarization of the tissue with AP propagation, in the case of fixed stretch. To support this point, in Figure 6 we provide a representative sequence of point-wise activations delivered on our simplified 2D domain and mimicking the experimental protocol conducted in Figure 1 for a selected parameter choice, i.e., (D1,Gs)=(-0.75·10-4,0). In this case, the AP excitation wave propagates differently according to the applied stretch state, both horizontal and vertical displacement and traction. In addition, the computed CVs change similarly to what observed in Figure 2. We remark that such a comparison with experimental observations is purely qualitative and does not represent a definitive validation of the model.

We further investigate the strength of MEF coupling effects. In particular, we want to determine which specific contribution (stretch-activated currents or stress-assisted diffusion) exhibits a better match against experimental evidence, and for this we assess changes in the S1-S2 stimulation protocol. In practice, in order to induce a spiral wave on an excitable tissue, one typically generates a planar electrical excitation (S1), followed by a second broken stimulus (S2) during the repolarization phase of the S1 wave, the so called vulnerable window (Karma, 2013). In our case, we selected a reduced set of MEF parameters (D₁, G_s) indicated in Table 3 as A,B,C,D. These values are motivated by the results from Figure 5. In particular, we select only the parameter combinations that produce either a unique decrement or increment of CV.

Figure 7 shows the different dynamics obtained via the S1-S2 protocol for the four different sets of MEF parameters. The first column is set at 100ms from the S1 stimulus for all the combinations, while the remaining frames are selected to highlight the elicited behavior. As a result, we observe that the deformation state of the tissue influences the overall dynamics differently. The first column highlights the variability in the AP wavelength, representing the spatial extension of the activation wave, which is due to the different repolarization states of the tissue induced by stress-assisted diffusion and stretch-activated currents. In particular, the AP wavelength varies as >6.2cm for case A, = 6.2cm for case B, and <2cm for cases C, D. In fact, when the G_s contribution is present, the excitation wave is much reduced with respect to the profiles generated with the electrophysiological three-variable model (1) and fine-tuned on experimental data. Such an effect is not present when G_s = 0.

Secondly, cases A and B (that is, where only D₁ is activated) provide a similar behavior for spiral onset and case B shows the expected reduction in CV. Contrariwise, cases C and D (where also the contribution of G_s is present) induce much more complex dynamics, not expected in an isotropic medium. In particular, case C leads to a wave break and multiple spirals generation at the S2 stimulus that eventually collide and result in a single spiral wave. On the other hand, case D shows a more stable behavior generated by the presence of D₁.

In addition, Table 3 also provides the minimum and maximum delay for the S2 stimulation (vulnerable window) allowing to induce a spiral wave in the uniaxially stretched tissue. It is evident that the presence of SAC reduces the minimum S2 stimulation time, tS2min, by about 100ms with respect to the other cases and slightly increase the overall time span of the vulnerable window. Such a variation is motivated on the additional reaction current induced by the presence of Isac(λ, u) everywhere in the medium, but it is not expected from the experimental isochrones provided in Figure 1.

To further corroborate this analysis, we provide in the top panels of Figure 8 an additional sequence referring to the combination (D1,Gs)=(-1.5·10-4,0.25) in the case with static displacement boundary conditions, which falls in the range where no CV wave was measured. As anticipated, an excessive contribution due to SAC elicits extra activations where the stretch is maximum, i.e., at the corners of the domain. This particular behavior is not obtained when the stress-assisted contribution D₁ is very high. Next, the bottom panels of Figure 8 show results using the combination (D1,Gs)=(-0.75·10-4,0.125), which allows the quantification of CV but can eventually lead to spiral breakup and non-sustainability of the arrhythmic patterns due to the mechanical state of the tissue (corresponding to the case of dynamic traction, described below). This is a representative example of the key importance of boundary conditions and how MEF effects could be effectively translated into clinical studies.

Finally, we turn to the analysis of meandering for the spiral tip for long run simulations (4s of physical time) comparing the four selected sets of parameters A,B,C,D in combination with static/dynamic–displacement/traction boundary conditions. In particular, we initiate the spiral wave via the S1-S2 stimulation protocol as discussed in the previous section, in absence of any mechanical loading such to start from the same initial conditions for each selected case. After spiral onset and stabilization (namely, for t > t₂ = 250ms), we apply the following four different loadings:

Static displacement: uniaxial displacement φ~=[0.1L,0]T applied on the right boundary while keeping the left one clamped (Figure 9A).

For each mechanical loading, panels in Figure 9 show the trajectories of the spiral tip for the four MEF parameters combinations. Two important aspects are worthy of attention.

First, for each combination of the mechanical loading, the presence of the stress-assisted conductivity D₁ tends to stabilize the meandering (see black and green traces). This behavior is particularly evident in Figure 9C where the combination D1=-0.75·10-4,Gs=0 results into a localized core, while the case D₁ = 0, G_s = 0 presents a circular, but slightly drifting core. Consequently, local stress-based heterogeneities appear in the medium when D₁ is different from zero, leading to pinning-like phenomena also observed in Cherry and Fenton (2008), Cherubini et al. (2012), Jiménez and Steinbock (2012), and Liu et al. (2013). Moreover, these conditions are associated with an ellipsoidal shape of the core underlying the effective anisotropy induced by the stress-assisted coupling. All these observations agree with the conclusions from the extended analysis conducted on the chosen AP model in the original work from Fenton and Karma (1998b).

Secondly, when also SAC is present, the spiral meandering is unpredictable and strongly dependent on the applied boundary conditions (see blue and red traces). In this scenario, it is interesting to note that static loading induces a simple meandering which eventually pushes the spiral wave out from the domain (see Figure 9C), whereas dynamic conditions dictate a chaotic behavior that makes the spiral either to explore the whole domain, or to exit it. These patterns seem to be extreme conditions of hyper-excitability not expected in a two-dimensional isotropic medium (Fenton and Karma, 1998a; Fenton et al., 2002).

Finally, we highlight the symmetry of the observed behavior according to the clockwise or counterclockwise rotation of the spiral. This particular analysis is provided in Figure 10 and further links the excitation dynamics to the mechanical features. The different traces refer to the spiral core meandering observed for a dynamic uniaxially stretched case with MEF parameters D₁ = 0, G_s = 0.125 and initiated via the S1-S2 stimulation protocol: case (a) compares a clockwise and counterclockwise spiral propagation; case (b) shows a counterclockwise spiral core initiated from the top (red) and bottom (blue) case. Corresponding sequences are also shown as side panels. This result is limited to the simplified nature of the domain adopted, i.e., 2D isotropic. A more realistic computational domain, embedding fiber directionality and tissue thickness, would show more involved dynamics in a complex spatiotemporal and clinical relevant perspective.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.