Competing Mechanisms of Stress-Assisted Diffusivity and Stretch-Activated Currents in Cardiac Electromechanics

We numerically investigate the role of mechanical stress in modifying the conductivity properties of cardiac tissue, and also assess the impact of these effects in the solutions generated by computational models for cardiac electromechanics. We follow the recent theoretical framework from Cherubini et al. (2017), proposed in the context of general reaction-diffusion-mechanics systems emerging from multiphysics continuum mechanics and finite elasticity. In the present study, the adapted models are compared against preliminary experimental data of pig right ventricle fluorescence optical mapping. These data contribute to the characterization of the observed inhomogeneity and anisotropy properties that result from mechanical deformation. Our novel approach simultaneously incorporates two mechanisms for mechano-electric feedback (MEF): stretch-activated currents (SAC) and stress-assisted diffusion (SAD); and we also identify their influence into the nonlinear spatiotemporal dynamics. It is found that (i) only specific combinations of the two MEF effects allow proper conduction velocity measurement; (ii) expected heterogeneities and anisotropies are obtained via the novel stress-assisted diffusion mechanisms; (iii) spiral wave meandering and drifting is highly mediated by the applied mechanical loading. We provide an analysis of the intrinsic structure of the nonlinear coupling mechanisms using computational tests conducted with finite element methods. In particular, we compare static and dynamic deformation regimes in the onset of cardiac arrhythmias and address other potential biomedical applications.


INTRODUCTION
Cardiac tissue is a complex multiscale medium constituted by highly interconnected units, cardiomyocytes, that conform a so-called syncitium with unique structural and functional properties (Pullan et al., 2005). Cardiomyocytes are excitable and deformable muscular cells that present themselves an additional multiscale architecture in which plasma membrane proteins and intracellular organelles all depend on the current mechanical state of the tissue (Salamhe and Dhein, 2013;Schönleitner et al., 2017). Dedicated proteic structures, such as ion channels or gap junctions, rule the passage of charged particles throughout the cell as well as between different cells and they are usually described mathematically through multiple reaction-diffusion (RD) systems (Cabo, 2014;Dhein et al., 2014;Kleber and Saffitz, 2014). All these coupled nonlinear and stochastic dynamics, emerge then to conform the coordinated contraction and pumping of the heart Land and et. al., 2016;Quarteroni et al., 2017). During the overall cycle, the mechanical deformation undoubtedly affects the electrical impulses that modulate muscle contraction, also modifying the properties of the substrate where the electrical wave propagates. These multiscale interactions have commonly been referred in the literature as the mechano-electric feedback (MEF) (Ravelli, 2003). Experimental, theoretical and clinical studies have been contributing to the systematic investigation of MEF effects, already for over a century; however, several open questions still remain (Quinn et al., 2014;Quinn and Kohl, 2016;Land et al., 2017;Sack et al., 2018). For example, and focusing on the cellular level, it is still now not completely understood what is the effective contribution of stretch-activated ion channels and which is the most appropriate way to describe them. In addition, and focusing on the organ scale, the clinical relevance of MEF in patients with heart diseases remains an open issue (Orini et al., 2017), and more specifically, how MEF mechanisms translate into ECGs (Meijborg et al., 2017) and what is the specific role of mechanics during cardiac arrhythmias (Christoph et al., 2018).
The theoretical and computational modeling of cardiac electromechanics has been used to investigate some key aspects of general excitation-contraction mechanisms. For instance, the transition from cardiac arrhythmias to chaotic behavior, including the onset, drift and breakup of spiral/scroll waves (Panfilov and Keldermann, 2005;Bini et al., 2010;Keldermann et al., 2010;Dierckx et al., 2015), pinning and unpinning phenomena due to anatomical obstacles (Cherubini et al., 2012;Hörning, 2012;Chen et al., 2014), as well as the multiscale and stochastic dynamics both at subcellular, cellular and tissue scale (Trayanova and Rice, 2011;Hurtado et al., 2016;Land et al., 2017). However, the formulation of MEF effects into mathematical models has been primarily focused on accounting for the additive superposition of an active and passive stress to stretch-activated currents (Panfilov and Keldermann, 2005). Recent contributions have advanced an energy-based framework for the comparison of active stress, stretch-activated currents and inertia effects (Cherubini et al., 2008;Ambrosi and Pezzuto, 2012;Rossi et al., 2014;Costabal et al., 2017). These works further highlight the role of mechanics into the resulting heart function at different temporal and spatial scales.
In order to further motivate our theoretical developments, we provide an experimental representative example of the strong MEF coupling in cardiac tissue, observable on the macroscale. The data shown in Figure 1 were obtained via dedicated fluorescence optical mapping applied on a pig right ventricle (the experimental procedure has been previously described in Fenton et al., 2009;Gizzi et al., 2013;Uzelac et al., 2017).
After motion suppression via blebbistatin, the perfused tissue was electrically stimulated via an external bipolar stimulator with strength twice diastolic threshold. An excitation pulse with constant pacing cycle length of 1 s was delivered within the field of view (red spot in Figure 1) for several seconds (reaching a steady-state configuration) and for three different mechanical loading conditions on the same wedge: (a) free edges, (b) static uniaxial horizontal stretch, (c) static uniaxial vertical stretch with respect to a prescribed tissue orientation. The figure displays the underlying structure with clear evidence of the deformed tissue architecture, isochrones of electrical activation for a representative stimulus, and a sequence of spatial activation maps, where the colors indicate the level of activation-Action Potential (AP). Since in this proof of concept setup active contraction is inhibited by blebbistatin, these experiments clearly indicate that an additional degree of heterogeneity and anisotropy appears in the tissue and affects the AP excitation wave due to the intensity and direction of the externally applied deformation. In addition, this behavior does not correspond to a mere linear mapping from the reference to the deformed configuration (as a visual scaling of the image would easily show), but one observes that mechanical deformations induce higher, nonlinear and non-trivial anisotropies and heterogeneities in the tissue.
To better characterize such features, in Figure 2 we provide an extended analysis of the local conduction velocity (CV) thorough histogram plots measured as follows: • we identify wavefront isochrones at 50% of depolarization for eleven consecutive frames at 2 ms each (this produces ten consecutive measures of CV per direction selected); • we compute the contour normal direction and the corresponding distance between consecutive isochrones; • we measure the local CV for all the computed normal directions, along the isochrone path and for seven consecutive action potential activations at constant pacing cycle length of 1 s; • we exclude the extreme values from the histogram to take out spurious results, e.g., boundary effects.
The chosen methodology allows to represent tissue heterogeneity, provides a robust measure of the local CV distribution characterizing the underlying ventricular structure, and homogenizes physiological beat-to-beat variabilities. We summarize the results of such an extended analysis in Table 1, distinguishing between the three loading cases as described in Figure 1, providing sample size and statistical features of the computed CV histogram distribution, i.e., mean and median. We also provide the box plot representation of the obtained distributions for the three stretch states, respectively, to further highlight dispersion of the measured velocities. Every single feature in the study confirms a slower conduction velocity under stretch, and this behavior is full agreement with previous studies (Ravelli, 2003). Also, in Figure 3 we demonstrate that the tissue is at steadystate for the selected stimulation rate providing a quantitative comparison of the spatial and temporal activation sequences. In particular, after several activations (> 5), beat n and beat FIGURE 1 | MEF observed in pig right ventricle via fluorescence optical mapping. From top to bottom, we provide: underlying tissue structure in reference (A) and deformed (B,C) states; activation isochrones each 4 ms originating from the stimulation point (red spot in the field of view-the bar indicates a length of 1 cm), and activation sequences. The three cases refer to no-stretch (A), static horizontally (B), and vertical (C) stretch in the directions indicated by the yellow arrows. The sequence of spatial activation uses the color code scaled to the AP level (yellow/green-high/low). Selected frames highlight the anisotropy induced by stretch. The outer black region is the noisy area not useful for the field of view.
FIGURE 2 | CV histograms measured on tissue wedges for three different loading states overlapping local measures for seven consecutive activations at constant pacing cycle length of 1 s. All the normal directions to the AP propagation are considered as indicated by orange arrows on a representative isochrone contour. The box plot of the distribution is provided as inset for the three histogram, respectively, highlighting the amount of dispersion and the reduction of CV under stretch (see Table 1 for details). Cut-off of spurious values is set at 0.05 and 1.3 m/s. n + 10 are shown for a selected frame in terms of normalized AP distribution and its spatial difference, as well as comparing the time course of two consecutive activations (B1, B2) for a representative pixel under the field of view. In both cases, the spatio-temporal differences recorded are within the physiological variability of a ventricular wedge, and the tissue shows a steadystate regime which is considered at resting state for the numerical model. Clear MEF effects evidenced in the previous experimental exercise suggest the incorporation of deformation and stress into the conduction properties of the cardiac tissue itself. The preliminary character of the proposed minimal model implies that we do not take into account the intrinsic structural variability of the tissue, but we stress that these effects will be investigated in future validation works. Accordingly, as a base line model, in the present study we will adapt the formulation recently proposed in Cherubini et al. (2017) and designed for general purpose stress-diffusion couplings. Doing so will allow us to readily and selectively incorporate two main MEF-related mechanisms into the computational modeling of cardiac electromechanics: (i) stretch-activated currents (SAC) and (ii) stress-assisted diffusion (SAD). The first paradigm relates the deformed mechanical state to the excitability of the medium via additional reaction functions (ionic-like currents); whereas the second one collects the homogenized effects of the deformation field on the diffusion processes originating the spatio-temporal patterns of the membrane voltage.
Within such a framework, we expect stretch-activated currents and stress-assisted diffusion to counterbalance each other by locally enhancing tissue excitability as well as smoothing the excitation wave according to the mechanical state of the tissue. In particular, since an external loading activates SAC at locations where the stretch is high and, at the same time, induces an heterogeneous and anisotropic diffusion tensor via the SAD mechanisms, our study focuses on the role of different mechanical boundary conditions in affecting action potential propagation and onset of arrhythmias. Accordingly, these two MEF mechanisms will be studied numerically in terms of three basic lines. First, by conducting a parametric analysis of the competing nonlinearities such to identify the limits of applicability of the proposed models. In particular, we select in the SAD mechanisms the most reliable modeling approach able to reproduce the experienced conduction velocity reduction upon an applied static loading state. Then, by performing a selective investigation of spiral onset protocols we will characterize the additional nonlinearities that arise due to MEF. Here we identify the different time span of the vulnerable window obtained via an S1S2 excitation protocol. Finally, by means of long-run analyses of arrhythmic scenarios, we compare and contrast static and dynamic displacement and traction loadings  Figure 1]. The first two rows show the spatial distribution of the normalized voltage for beat n and beat n + 10 with the corresponding difference in the third row (color code is indicated). The last row indicates the time course of a representative pixel in the center of the field of view for two consecutive beats n and n + 10 with the corresponding difference provided in the red trace. on a two-dimensional, idealized tissue slab. In this regard, we show how spiral core meandering results highly affected by the mechanical state and becomes unstable when SAC and SAD parameters are stronger.
Our results highlight several interesting conclusions regarding the propagation of the excitation wave in the presence of two competing MEF effects. These findings call for novel and additional experimental investigations. Finally, we provide a thorough discussion of the applicability of the proposed modeling approach and its extensions toward more realistic and multiphysics scenarios.

METHODS
The classical stress-assisted formulation proposed in Aifantis (1980) was developed in the context of dilute solutes in a solid. A similarity exists between this fundamental process and the propagation of membrane voltage within cardiac tissue. Indeed, on a macroscopically rigid matrix, the propagating membrane voltage can be regarded as a continuum field undergoing slow diffusion. Here we consider a similar approach (developed in Cherubini et al., 2017) which generalizes Fick's diffusion by using the classical Euler's axioms of continuously distributed matter.
In particular, the balance of momentum can be imposed such to ensure frame invariance, a property of high importance in mechanical applications (Tadmor et al., 2012). We also assume quasi-static conditions for the continuum body, such that its macroscopic response is, in principle, independent from the diffusion process. On the contrary, the diffusion process will strongly depend on the mechanical state of the tissue.

Continuum Electromechanical Model
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 = det F iI the Jacobian of the map and with C IJ = F kI F kJ and B ij = F iK F 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 iJ a 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 1 and I 2 rule the deviatoric response of the medium, the third invariant I 3 quantifies volumetric changes of the material, while the fourth pseudo-invariant I 4 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 stressassisted 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 (2a), i.e., [d ij ∂u/∂x j ]n i = 0, where n i is the outward normal on the domain boundary. System (2) describes the propagation of a normalized dimensionless membrane potential, which can be mapped to physical quantities as u = 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 (2a), 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 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 I max ext . 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 ρ 0 , ρ and dV, dv are the densities and volumes of the solid in the undeformed and deformed configurations, respectively.  Fenton and Karma (1998b) and Cherubini et al. (2017).
Time units are ms, length is cm, the termḡ fi is in mS/cm 2 , dimensional voltages are in mV, and stiffness in MPa. Square brackets indicate range of parameter variability, and the rightmost column specifies initial conditions for a resting tissue.
In Equation (3b),φ ϕ ϕ(t) is a known (possibly time-dependent) displacement and in Equation (3c),t i (t) is a (possibly timedependent) 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 quasistatic 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 1 and c 2 ; 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 (4b), where the constant k Ta modulates the amplitude of the active stress contribution, while ǫ(u) is a contraction switch function: ǫ(u) = ǫ 0 if u < 0.005, and ǫ(u) = 10ǫ 0 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 0 represents the usual diffusion coefficient for isotropic media, i.e., diffusivity = [L 2 T −1 ], while D 1 and D 2 introduce the impact of mechanical stress through linear and nonlinear contributions, respectively, on the diffusive flux. Accordingly, D 1 and D 2 have units of [L 2 T −1 P −1 ] and [L 2 T −1 P −2 ], respectively. We also remark that Equation (4c) reduces to the characterization of the classical diffusion equation for D 1 ≡ D 2 = 0.
Finally, Equation (4d) describes the stretch-activated current contribution (which is usually adopted as the sole MEF effect). The term I sac (λ, 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, u sac 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 tip X I + y tip Y 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 −4 .

Numerical Approximation
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 T h of (0) into triangles T of diameter h T , where h = max{h T : T ∈ T h } is the meshsize. The finite element spaces mentioned above are defined as (see e.g., Quarteroni and Valli, 1994) H h : = {ψ ∈ H 1 ( (0)) : ψ| T ∈ [P 2 (T)] 2 ∀T ∈ T h , and ψ = 0 on Ŵ D (0)}, Q h : = {q ∈ L 2 ( (0)) : q| T ∈ P 1 (T) ∀T ∈ T h }, for the case of clamped boundaries at Ŵ D (0). Let us also construct an equispaced partition of the time domain 0 = t 0 < t 1 = t < · · · < t M = t max . The coupled problem is solved sequentially between the mechanical and electrochemical blocks. A description of the needed computations at each time step t n is as follows: Step 1: From the known values u n h , v n h , w n h , T n a,h , D n h , λ n h , find This scheme for the electric/activation system is given in a first-order semi-implicit form: the nonlinear reaction terms and the coupling stressassisted 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 T n+1 a,h computed in Step 1 of this iteration, solve the nonlinear elasticity equations 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 −6 . 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 −5 . The residual computation for the mechanical problem also contains the terms arising from timedependent 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.2L sin 2 (π/400 t), 0 T .
All tests are conducted using a two-dimensional slab of dimensions L × L = 6.2 × 6.2 cm 2 , 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.1 ms, h = 0.062 cm, 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 (Alnaes et al., 2015).

RESULTS
In the following, we adopt a parametric setup fitted for the modified Beeler-Reuter model (Equation 2), while selectively changing MEF parameters (D 1 , G s ). This choice provides a reference, unloaded, model configuration with constant CV of 0.42 m/s and a circular meandering for a free spiral on a homogeneous and isotropic domain. Such values deviate as the MEF coupling is activated.

Conduction Velocity Analysis
We start analyzing the parameter space associated to the two MEF contributions in our model. That is, the stress-assisted coefficients D 1 , D 2 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 2 to be one order of magnitude smaller than D 1 . 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 1 , 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 1 = 0, G s = 0), we observe a net increase of CV for (D 1 = 0, G s > 0) while we recover CV decrements for (D 1 < 0, G s = 0). This specific FIGURE 4 | Example of structured mesh employed in the computational results. The grid is displayed on the deformed configuration when the domain is subject to traction (arrows) and fixed displacement (lines) boundary conditions, and a zoom exemplifies the mesh resolution for a rather coarse spiral front. 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 1 . Second, for the particular combination (D 1 < −10 −4 , G s > 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., (D 1 , G s ) = (−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.

S1-S2 Excitation Protocol
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  Table 3) are highlighted together with two additional cases in which CV was not recorded. On the right, three consecutive time frames of the activation are selected. 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 1 , 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 100 ms 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.2 cm for case A, = 6.2 cm for case B, and < 2 cm for cases C, D.
FIGURE 7 | S1-S2 stimulation protocol applied on a static uniaxial stretched configuration for different combinations of MEF parameters (D 1 , G s ) as provided in Table 3. The color code refers to normalized dimensionless membrane potential, u, (blue-red mapped to [0-1]). Selected time frames are provided in the subpanels.
Frontiers in Physiology | www.frontiersin.org 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 (2) 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 1 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 1 .
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, t min S 2 , by about 100 ms 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 I sac (λ, 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 (D 1 , G s ) = (−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 1 is very high. Next, the bottom panels of Figure 8 show results using the combination (D 1 , G s ) = (−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.

Spiral Drift and Effects due to Boundary Conditions
Finally, we turn to the analysis of meandering for the spiral tip for long run simulations (4 s 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 2 = 250 ms), 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). • Dynamic displacement: uniaxial time-dependent displacementφ ϕ ϕ(t) = 0.1L sin 2 (π/400 t), 0 T applied on the right boundary while keeping the left one clamped ( Figure 9B). • Dynamic traction: uniaxial time-dependent forcet i (t) = t max sin 2 (π/400 t) applied on the left and right boundaries while keeping the bottom side clamped ( Figure 9D).
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 1 tends to stabilize the meandering (see black and green traces). This behavior is particularly evident in Figure 9C where the combination D 1 = −0.75 · 10 −4 , G s = 0 results into a localized core, while the case D 1 = 0, G s = 0 presents a circular, but slightly drifting core. Consequently, local stress-based heterogeneities appear in the medium when D 1 is different from zero, leading to FIGURE 9 | Tip trajectories for four combinations of MEF parameters (D 1 , G s ) (see Table 3 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 twodimensional 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 1 = 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.

CONCLUSION
We have advanced a minimal model for the electromechanics of cardiac tissue, where the mechano-electrical feedback is incorporated through two competing mechanisms: the stretchactivated currents commonly found in the literature, and the stress-assisted diffusion (or stress-assisted conductivity) recently proposed by Cherubini et al. (2017). Both the electrophysiology and the mechanical response adopt a phenomenological simplified description, but a preliminary validation is provided through a set of numerical simulations that agree qualitatively with a set of experimental data for pig right ventricle.
The implications of the intensity and degree of nonlinearity assumed for the stress-assisted diffusion effect are studied from the viewpoint of changes in the conduction velocity and the dynamics of spiral waves in simplified 2D domains. Multiple electrical stimulations protocols and non-trivial mechanical loadings have been investigated highlighting the strong coupling due to the different MEF contributions. The analysis supports the hypothesis that the simplistic formulation adopted for stretch-activated currents seems to deviate from the experimental evidence, in line with recent contributions addressing the coupled modeling of SACs and stretch-induced myofilament calcium release at the myocyte level (Timmermann et al., 2017). On the other hand, in a homogenized setting, the stress-assisted diffusion formulation produces a series of interesting phenomena that qualitatively match heterogeneities and anisotropies observed during mechanical stretching of pig right ventricle via fluorescence optical mapping.
Limitations of the present work are partially linked to the phenomenological approach adopted to describe the complex multiscale mechanisms intrinsic in the cardiac tissue and partially due to the simplified computational domain. In this regards, we aim at investigating more reliable stretch-activated current formulations leading to alternans behaviors (Galice et al., 2016) within a multiscale mechanobiology perspective (Nava et al., 2016;Stålhand et al., 2016;Cyron and Humphrey, 2017) and tacking into account the intracellular calcium cycling influenced by mechanical stretch, because all these effects have been proposed as concurring mechanisms of arrhythmogenesis within the heart. From the mechanical point of view, we mention as main limitation the adoption of a simplified isotropic hyperelastic material model which can be generalized to more complex and reliable formulations. This will include, for example, active strain anisotropies, muscular and collagen fiber distributions in an orthotropic mechanical framework that the authors have been extensively developing during the last decade (Cherubini et al., 2008;Nobile et al., 2012;Gizzi et al., 2015Gizzi et al., , 2016Gizzi et al., , 2018Pandolfi et al., 2016). Such a generalization will maintain the nature of the present theoretical framework in terms of MEF competing effects. In this line, we also aim to generalize our theoretical and computational approach toward intrinsic multiscale and multiphysics mechano-transduction problems (Weinberg et al., 2017;Lenarda et al., 2018), e.g., the uterine smooth muscle activity (Young, 2016;Yochum et al., 2017) or the intestine biomechanics activity (Pandolfi et al., 2017;Brandstaeter et al., 2018) by implying the usage of network approaches (Giuliani et al., 2014;Robson et al., 2018) and data assimilation procedures (Barone et al., 2017). In addition, the investigation of the complex spatiotemporal dynamics, chaos control and multiphysics couplings in excitable systems (see e.g., Hörning et al., 2017;Christoph et al., 2018) can be emphasized within the proposed electromechanical framework by using realistic three-dimensional cardiac structures (Lafortune et al., 2012). We also mention implications of the proposed models in the mathematical study of general stress-assisted diffusion problems, as recently carried out in Gatica et al. (2018). Finally, we hope that the present contribution may open new experimental studies to translate the complex MEF phenomena into the clinical practice (Meijborg et al., 2017;Orini et al., 2017) identifying novel risk indices for cardiac arrhythmias .

ETHICS STATEMENT
All experiments conform to the current Guide for Care and Use of Laboratory Animals published by the National Institutes of Health (NIH Publication No. 85-23, revised 1996), and approved by the Office of Research and Integrity Assurance at Georgia Tech.

AUTHOR CONTRIBUTIONS
AL, AG, RR-B, CC, FF, and SF design of the study; AL, AG, and RR-B numerical methods and computational tests; FF experimental measurements; AL and AG Statistical analysis; AL, AG, RR-B, CC, FF, and SF Manuscript writing.

ACKNOWLEDGMENTS
This work has been supported by the Italian National Group of Mathematical Physics GNFM-INdAM; by the International Center for Relativistic Astrophysics Network ICRANet; by the London Mathematical Society through its Grant Scheme 4; and by the EPSRC through the Research Grant EP/R00207X/1.