Spatiotemporal Organization of Electromechanical Phase Singularities during High-Frequency Cardiac Arrhythmias

Ventricular fibrillation (VF) is a lifethreatening electromechanical dysfunction of the heart associated with complex spatiotemporal dynamics of electrical excitation and mechanical contraction of the heart muscle. It has been hypothesized that VF is driven by three-dimensional (3D) rotating electrical scroll waves, which can be characterized by filament-like electrical phase singularities (EPS). Recently, it was shown that rotating excitation waves during VF are associated with rotating waves of mechanical deformation. 3D mechanical scroll waves and mechanical filaments describing their rotational core were observed in the ventricles by using high-resolution ultrasound. The findings suggest that the spatiotemporal organization of cardiac fibrillation may be assessed from waves of mechanical deformation. However, the complex relationship between excitation and mechanical waves during VF is currently not understood. Here, we study the fundamental nature of mechanical phase singularities (MPS), their spatiotemporal organization and relation with EPS. We demonstrate the existence of two fundamental types of MPS:"paired singularities", which are co-localized with EPS, and"unpaired singularities", which can form independently. We show that the unpaired singularities emerge due to the anisotropy of the active force field, generated by fiber anisotropy in cardiac tissue, and the non-locality of elastic interactions, which jointly induce strong spatiotemporal inhomogeneities in the strain fields. The inhomogeneities lead to the breakup of deformation waves and create MPS, even in the absence of EPS, which are typically associated with excitation wave break. We exploit these insights to develop an approach to discriminate paired and unpaired MPS. Our findings provide a fundamental understanding of the complex spatiotemporal organization of electromechanical waves in the heart.


I. INTRODUCTION
The beating of the heart is initiated by nonlinear waves of electrical excitation, which propagate through the cardiac muscle and trigger a release of intracellular calcium, which in turn triggers contractions of cardiac muscle cells.During severe heart rhythm disorders, such as atrial or ventricular fibrillation, the electrical excitation degenerates into multiple disorganized, asynchronous electrical waves leading to irregular cardiac muscle contractions [1][2][3][4].Both spiral-shaped re-entrant waves and focal waves are thought to underlie those phenomena during highfrequency cardiac arrhythmias.To date, however, the three-dimensional (3D) electrical wave phenomena that evolve rapidly within the thickness of the heart muscle have never been imaged in full.While fluorescence imaging provides high-resolution measurements of vortex-like rotating spiral waves on the heart surface, the underlying 3D dynamics remain elusive, and, based on numerical simulations, are conjectured to take on the shapes of scroll vortex waves or 3D spiral waves [5][6][7][8][9][10][11], which are generic self-organized nonlinear wave structures of excitable media [1-4, 12, 13].
In a recent study, it was shown that the dynamical processes underlying ventricular fibrillation can be characterized by coupled electrical and mechanical phase singularity dynamics [14].By using tri-modal voltage-and calcium-sensitive fluorescence imaging, high-speed 4D ul-trasound, and numerical motion analysis, it was shown that the rapidly contracting, fibrillating heart muscle exhibits vortex-like rotating strain fields, which are produced by electrical action potential and calcium spiral vortex waves.Consequently, it was demonstrated on the heart surface that electrical phase singularities, which describe the cores of rotating action potential and calcium waves, appear in the vicinity of the cores of the rotating deformation patterns, which can equivalently be characterized by mechanical phase singularities.It was furthermore shown that electrical and mechanical phase singularities co-localize and exhibit similar dynamics in terms of numbers, trajectories, and life-times, providing for the first time evidence for the existence of electromechanical rotors.The data suggests that electrical and mechanical phase singularities are both produced during ventricular fibrillation by electromechanical scroll wave chaos, and that a better understanding of the nature of these electromechanical phase singularities could provide important new insights into the 3D spatiotemporal organization of cardiac fibrillation since electrical and mechanical phase singularities are closely related to each other.Lastly, it was shown using high-resolution 4D ultrasound, that a 3D filament-like structure of mechanical phase singularities evolves throughout the heart wall during fibrillation, suggesting that the spatiotemporal organization of electrical scroll wave vortex filaments may be inferred by the dynamics of mechanical vortex filaments inside the heart muscle.The findings open the path to use mechanical phase singularities to enhance the understanding of the mechanisms underlying cardiac fibrillation.However, despite the experimental advances, a fundamental understanding of their relationship remains largely lacking.
Naively, one would expect mechanical waves to be slaved to electrical excitation waves that cause contraction via the well-established excitation-contraction coupling mechanism, whereby calcium entry into the cell through L-type calcium channels following electrical excitation triggers calcium release from intracellular stores, which in turn activates the contractile machinery of the cell thereby generating an active force along the long axis of cardiomyocytes [15].Consequently, the spatiotemporal wave pattern of this active force is expected to generally follow closely the wave pattern of electrical excitation, except under extreme conditions where excitationcontraction coupling fails (e.g., , when calcium-inducedcalcium-release fails to keep up with electrical excitation at very high frequencies).When normal excitationcontraction coupling remains operative, excitation wave propagation can itself be influenced by mechanical contraction due to mechanoelectrical feedback [16][17][18].However, this feedback, which has been neglected in most modeling studies of cardiac arrhythmia mechanisms [1][2][3], does not cause the wave pattern of active force to deviate from the wave pattern of electrical excitation.In contrast, the resulting tissue deformation pattern, which is imaged by high-resolution ultrasound, may differ markedly from the active force pattern due to the non-locality of long-range elastic interactions, which can induce strain in regions of the heart muscle that are not actively contracting or, at the opposite, cause regions under large active force to be under small strain.This raises the basic question of whether mechanical phase singularities, i.e., spatiotemporal phase singularities of non-local strain fields, track electrical wave singularities, or exhibit a more complex spatiotemporal organization.
In this article, we address this question by investigating the strain field patterns generated by re-entrant and focal excitation-contraction waves in 2D and 3D cardiac tissue in computer simulations.Our main finding is that mechanical phase singularities, which we extract from phase maps obtained of strain fields via the Hilbert transform, can both co-localize with, or form away from, phase singularities of excitation waves.The former "paired" mechanical singularities only exist in the presence of reentrant waves while, remarkably, focal sources of waves or target patterns suffice to create the latter "unpaired" mechanical singularities, which can, therefore, exist both in the presence and absence of re-entrant waves.Further, we find that paired and unpaired singularities form by different mechanisms.Paired ones originate from phase singularities associated with the rotation of strain fields that remain essentially slaved to excitation waves near the core of electrical vortices.In contrast, unpaired ones originate from the breakup of deformation waves caused by the combination of long-range elastic interactions and the anisotropy of the active force field.
We exploit these insights to develop an approach to discriminate paired and unpaired mechanical phase singularities based on observations of their surrounding strain field.Our findings provide a fundamental understanding of the complex spatiotemporal organization of electromechanical wave activity during cardiac arrhythmias, and a theoretical basis for the interpretation of threedimensional imaging data of tachyarrhythmias obtained with high-resolution ultrasound for diagnostic and therapeutic applications.

II. METHODS
To explore the relationship between electrical and mechanical singularities, we use two different approaches.First, in section II A, we analytically construct simple 2D excitation wave patterns and compute the resulting strain fields (deformation waves) by numerically solving the equations of linear elasticity using the finite element method (FEM).We should note that there is a large body of literature dedicated to creating accurate representation of the coupled electromechanical response of the heart [19][20][21][22][23][24].However, since our goal is to elucidate the mechanism of mechanical singularity formation we reduce the problem to its most essential aspects.For the main presentation in this article, we use quasi-static for-mulation with homogeneous isotropic infinitesimal elasticity, with anisotropic active strains and no electromechanical feedback.We also limit our presentation in the main text to Poisson's ratio of ν = 0.4 which produces ∼ 10% volumetric change (see Fig. 3).However, to show the robustness of our findings, we provide an extended set of simulations using ν = 0.49, neo-hookean material, and transversely isotropic elasticity in the Supplementary Videos.Under the assumption that the local active force tracks the local electrical excitation, which is generally satisfied except when calcium release cannot keep up with very high-frequency electrical waves, we do not need to specify separately the electrical and mechanical signals.We prescribe directly the spatiotemporal active strain fields corresponding to standard excitation 2D wave patterns.As generic examples of reentrant and focal excitations, we consider both a rigidly rotating (i.e., non-meandering) Archimedean spiral wave and a target pattern generated by a localized time-periodic focal source of waves, respectively.This approach has the advantage that it provides the simplest possible setting to investigate the relationship of electrical and mechanical singularities in a perfectly elastic, nearly incompressible, material where the strain field is assumed to relax instantaneously to the active force field.In this approach, "electrical singularities" are simply the singularities of the prescribed active force fields, which are present and absent for spiral and target waves, respectively, and "mechanical singularities" are the singularities of the numerically obtained linear elastic strain fields generated by the active force fields.
Second, in section II B, we simulate electromechanical wave activity using a simplified two-variable ionic model of electrical excitation [25,26] and a simple phenomenological relaxational kinetic equation to relate the active force to electrical excitation [17].In addition, we compute the resulting strain fields by extending a mass-spring model (MSM) of standard elasticity [27] to incorporate the active force and describe a nearly incompressible viscoelastic material with inertia.This approach has the advantage of allowing us to efficiently explore the relationship of excitation-contraction and deformation waves in 3D anisotropic tissue in a parameter limit of the MSM where viscous and inertial effects are small and continuum elastic properties are similar to those obtained by quasistatic FEM solutions, with known continuum elastic properties that can be derived analytically from the MSM.
A. Two-dimensional FEM computations with analytically prescribed spiral and target wave active force fields We perform the computations using quasistatic linear elasticity, where we implemented active contraction by an analogy to thermal stress with the incorporation of (stress-free) eigen strains.Our use of active strains ensures the ellipticity of the elastic energy [28].For domain Ω ⊂ R 2 we write the elastic energy as: where C ijkl = λδ ij δ kl + µ(δ ik δ jl + δ il δ jk ) is the fourth order isotropic elastic constitutive tensor where δ ij is the Kronecker delta and for plane-stress λ = Eν/(1 − ν 2 ) and µ = E/(2(1 + ν)) are Lamé parameters, ij = (u i,j + u j,i )/2 is the linearized strain tensor (i.e., symmetric displacement gradient), 0 ≤ T a ≤ 0.15 is the imposed spatiotemporally varying active contraction field normalized by the elastic modulus E, and β is a tensor whose form given below depends on whether contraction is isotropic, as in the case of randomly oriented cardiomyocytes in a tissue culture, or aniostropic as in the case of heart tissue with aligned fibers.To obtain the quasistatic response, we find the displacement field u * that minimizes the elastic energy, Eq. ( 1).We assume cardiac tissue to be elastically isotropic in all cases studied here with i.e., Young's modulus E and Poisson's ratio ν = 0.4 to account for nearincompressibility of the tissue, which results in ∼ 10% maximum volume change in our FE simulations.This value of the Poisson ratio was selected to correspond to the value achieved by the MSM model given our choice of parameters (see Appendix D).To assess the effect of volume change on our findings, we duplicated our simulations for higher Poisson's ratio ν = 0.49 (See Supplementay Movies).We investigate two different 2D cases: Isotropic case.
The "isotropic" case mimics a thin quasi-2D tissue culture of randomly oriented cardiomyocytes with no preferred direction of conduction or active contraction.In this case, both the conduction and active strains are assumed to be isotropic.Furthermore, we set to produce an isotropic contraction of unit magnitude.Anisotropic case.
The "anisotropic" case represents normal 2D cardiac tissue with cardiomyocytes aligned along a common fiber axis where we assume the conduction along the fibers to be 5 times faster than conduction perpendicular to them.We further assume that the active strain only acts along the fiber.In this case, we choose the expansion tensor such that for a 1D fiber, the strains along the fiber converge to unit magnitude without any change of fiber cross-section.For a domain with fibers along the first coordinate axis, and assuming isotropic elasticity, this is achieved by choosing We note that the above equations ( 2) and (3) are special cases of a domain with arbitrary fiber orientation defined by the unit vector f with components (f 1 , f 2 ) with β defined by To perform the simulation, we explicitly impose the active contraction field T a using closed-form expressions in the form of spatially-diffuse Archimedean spirals and expanding ellipsoidal pulses (see Appendix A) in an L × L square domain discretized using 320 × 320 Q1 elements.We impose traction-free boundary conditions on the edges of the domain (x = ±L and y = ±L) and remove the rigid body modes (null-space of elasticity) from the obtained discrete set of equations.For our FEM numerical implementation, we use libMesh [29] for finite element book-keeping and PETSc [30,31] for linear algebra.
B. Three-dimensional simulations using a two-variable ionic model with active force generation coupled to a mass-spring model

Reaction diffusion model for cardiac excitation
In 3D simulations, we assume a slab of ventricular muscle with a defined fiber architecture.A two-variable reaction-diffusion model presented in [25,26] has been utilized to simulate the cardiac excitation.The magnitude of the active force is coupled phenomenologically to the voltage field to mimic the typical rise and fall of this force during an action potential [17].The equations that describe this model are: subjected to the Neumann boundary condition: In these equations, U and v are the dimensionless representation of the transmembrane voltage and slow current gate variables, respectively, τ U /τ v controls the relative abruptness of the excitation, D is the diffusivity tensor, and n is the unit normal vector to the boundary.In Eq. ( 7) (similar to Eq. ( 1)), T a represents the ratio of the active contraction field and elastic modulus.K t controls the maximum amplitude of the active contraction by coupling it to the trans-membrane voltage U , and χ(U ) is a step function that sets the time scale of the contraction period with respect to U .Details about functions f , g, and the relation between diffusivity tensor and the fiber architecture can be found in [26,32].For convenience, these relations and more details about Eq. 7 are briefly presented in Appendix B. Finally, with no loss of generality, we assume that the conduction is transversely isotropic, i.e., it is isotropic in the plane perpendicular to fibers and is faster along the fiber direction.We can therefore write D ⊥1 = D ⊥2 = D ⊥ where D ⊥1 and D ⊥2 are the diffusivity perpendicular to the fiber axis in each plane and D is the diffusivity along the fiber axis.This will simplify the local electric current to where f is the unit vector that is locally parallel to the fiber.Here we consider the anisotropic cases without fiber rotation where f is spatially uniform and with fiber rotation where f rotates along one axis as in [26,32].The set of parameters that are used in this article is listed in Table I.

Three-dimensional lattice mass-spring model
In the 3D simulations, the elasticity of the tissue was modeled by extending the mass-spring model (MSM) of Refs.[27] to cardiac tissue.Our use of the discrete MSM model (as opposed to discretizing a set of continuum-level equations), allows us to easily integrate the equations of motion on massively parallel GPUs.The schematic of a single unit cell of the cubic lattice with edge length a is shown in Fig. 1 and the tissue is constructed by stacking several unit cells.The extension to cardiac tissue was achieved by 1) adding Kelvin-type dampers (Fig. 1(b)) to the MSM to account for the viscoelastic behavior of the tissue, which produces damping forces proportional to the relative velocity between two masses, 2) a volumetric penalty force to shift the elastic properties towards the incompressible limit (Fig. 1(c)), and 3) by adding the active contractile force (Fig. 1(d)).For small and slowly spatially varying deformations (i.e., varying on a scale much larger than a) and in the absence of the volumetric penalty forces, the MSM reduces in the continuum limit to standard isotropic linear elasticity with Lamè constants λ = µ = k/a characteristic of a compressible material [27].
However, similar to other biological tissues, the myocardium contains mostly water and can be considered as a nearly incompressible material [19].To account for this near-incompressibility in our model, we introduce a uniform penalty pressure p to keep the volume constant.This pressure is applied to all faces of a unit cell and transmits a force on each mass at the corners of the cell.A schematic of the volumetric forces for the case that the volume decreases is shown in Fig. 1(c).The magnitude of the volumetric penalty force is chosen as where V 0 = a 3 and V are the initial and current volume of each cubical element, respectively.To ensure the balance of the internal forces, we apply the volumetric forces / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = "      on each pair of diagonally opposite nodes along the line connecting them.This volume preservation method can be shown to modify the bulk stiffness of the MSM.In the presence of the volumetric pressure, p, the first Lamé constant is increased to λ = k/a + p while the shear modulus remains unchanged µ = k/a.Therefore we can write Young's modulus and Poisson's ratio of the system as: In the simulations, we chose p = 15k/(2a) corresponding to a Poisson's ratio ν = 17/38 ≈ 0.45 close to the incompressible limit.We should note that while the value of penalty forces (10) vary spatially, the underlying model remains isotropic and homogeneous.Next, we incorporate the active contraction in our model.The active tension can be incorporated in two ways: (i) using active strains by dynamically setting the rest length of passive strains or (ii) using active stresses.Due to its ease of use in the MSM model, we opt for the latter option.We should note that these two options are equivalent only at small strains and can produce quantitatively different results at finite strains [28].We assume that contraction only occurs along the local fiber axis.In each unit cell, we assume that there is a pair of in-ternal active forces oriented in the direction of the fiber and choose the magnitude of these forces as |f a | = Ta Ea 2 where Ta = 8 i=1 T i th cell a /8 represents the average contractile field in the cell.The active force contribution from each cubical cell T i th cell a at a corner is calculated using bilinear interpolation based on the position where the fiber axis intersects the plane containing the mass, as described previously [33,34], and shown in Fig. 1(d).Finally, the equation of motion for the i th mass is given by where the various forces shown are in Fig. 1.They include the static spring forces f s = k×(relative displacements of masses linked by springs) and associated damping forces f d = c×(relative velocities) for the 18 dampers connected to each mass, and the active forces f a and volumetric forces f v defined above for the cubical elements connected to each mass.A detailed study of the mechanical properties of the MSM has been conducted and presented in Appendix D.
To study the electrical and mechanical singularities in heart tissue, we perform the simulations on a cubic lattice with 150 × 150 × 40 grid points.Adjusting for the material properties, the lattice corresponds to a 4.5 × 4.5 × 1.2 cm 3 slab with a unit cell spacing ∆x = ∆y = ∆z = a = 0.3 mm.To integrate the electrophysiology equations (Eqs.( 5)-( 7)) we use an explicit Euler method.To create a spiral wave initial condition, we use the standard two-stimulus protocol that consists of first creating a traveling plane wave along the y direction by exciting the tissue uniformly on the y = 0 plane, and then depolarizing half of the tissue (x < 2.25 cm) at t = 0.35 s.In addition, we integrate the equations of motion of the MSM using a standard velocity Verlet algorithm.
In a human heart with a density close to water ρ ≈ 1000 kgm −3 and Young's modulus ≈ 125 MPa [18], the elastic wave speed can be crudely estimated to be C ≈ 353 ms −1 .Furthermore, the period of excitation T spans the range 0.1 to 1.0s with the lower and upper bounds corresponding to an electrical rotor and a normal heartbeat, respectively.It follows that the distance CT traveled by elastic waves during a contraction period is much larger than the characteristic size L of the heart of a few centimeters.This scale separation (L CT ) makes it possible to increase the computational efficiency of the time integration algorithm by reducing the elastic wave speed.Therefore, instead of setting Young's modulus to its physical value, we set the characteristic elastic wave speed C = E/ρ equal to 5 ms −1 such that the inequality L CT still holds, where ρ = m/a is the tissue density.
The passive response of cardiac myocyte exhibits hysteresis cycles indicative of a viscoelastic behavior [35].While these effects can create small quantitative differ-ences compared to the non-viscoelastic models there is no evidence that they can significantly change the response of the tissue which is primarily driven by the active stresses [23].As such in our implementation we neglect such effects and we choose the damping constant such that the length of the system is much smaller than the spatial decay length of linear elastic waves, or L (CT )/(πδ), where tan(δ) = E l /E s with E l and E s are the loss and storage moduli, respectively, of our viscoelastic medium.In addition, we choose the elastic diffusion constant D e = η/ρ = 11 cm 2 s −1 where η = 5c/(2a) is the Kelvin-Voigt damping constant.We present a more detailed justification for these choices of parameters in Appendix C. The parameters for the MSM model are summarized in Table I.Fig. 1 implies that the sum of moments and forces in each cube and consequently in the full system is always zero.Therefore, all the forces in the system are internal and no additional mechanical boundary conditions need to be applied to simulate stress-free boundaries of the tissue.
We study two different fiber architectures: 1) all fibers being uniformly aligned in the x-direction and 2) fibers being organized in orthotropically stacked sheets, rotating about the z-axis through the thickness of the slab.For the second case, we assume that the fibers are parallel to the x-axis at the bottom (endocardium) and parallel to y-axis at the top surface (epicardium) with a total rotation angle of 90 • degrees.We first study spiral-shaped and circular-shaped focal electromechanical wave patterns in 2D isotropic tissues, see Figs. 2-6.Before we begin to present the result, we want to emphasize that the study conducted in this paper is not limited to a particular strain signal such as Tr( ) or xx or any particular direction, but it works for any signal that captures the mechanical response of the system.In this article we chose the strain along the fibers when there exists a fiber architecture (anisotropic case) or Tr( ) when there is not (isotropic case) to capture the mechanical response of the system.In the isotropic case, there is no particular fiber orientation (f i = 1 in Eq. ( 4)) and active tension occurs equally in all directions.Moreover, instead of modeling the electrical spiral by using reactiondiffusion kinetics, we directly derived the shape of the active contraction field T a analytically using a closed-form expression presented in Eq. ( 20) in Appendix A. Panels a) and b) in Fig. 2 show the contraction field T a (x, y) and the resulting deformation of the tissue displayed as strain (x, y) in material coordinates, respectively.Note that, because the contraction is isotropic, we chose the trace of the strain tensor Tr( )(x, y) = xx + yy to visualize the strain as a scalar-valued field, where negative values correspond to contracted tissue.In the isotropic case, it can generally be observed that the spiral shape emerging in the strain field matches the spiral shape in the active contraction field very well.The high correlation between the scalar-valued contraction field and a scalar-valued representation of the strain furthermore manifests when comparing time-series obtained from individual tissue segments: over several rotations of the spiral, each local rise and decline in active contraction T a (t) ∈ [0, 0.15] causes a corresponding shortening and elongation of the same tissue segment, which can be measured as a correspondingly oscillating strain signal s(t) = Tr( )(t) which aligns with the time-course of the contraction variable T a (t) for each and every material segment throughout the tissue.However, despite the strong correlation between contraction and resulting strain, the spiral shape in the strain field in (b) does not perfectly match the contraction spiral in (a), and scaling the fields would be insufficient to map one pattern onto another.Upon close inspection, the strain spiral pattern in Fig. 2(b) exhibits slight distortions and deviations, which can be observed best close to the tissue boundaries, c.f. Fig. 3(b).This behavior is to be expected because the overall strain field results from a superposition of several elastic phenomena which are all long-range effects interacting with each other over long distances throughout the entire medium: local contraction-induced deformations, such as contracted or elongated tissue regions pulling and pushing each other, and mechanical boundaries restricting the deformation.Note, that the strain fields are shown in Figs.2-6 are obtained at equilibrium and are fully relaxed.
Next, we use the Hilbert transform to assign a phase angle φ(t) and its complex amplitude A(t) to each deformation state of each material segment throughout the medium yielding a phase map φ(x, y, t) and an amplitude map A(x, y, t) as shown in panels (c) and (d) in Fig. 2, respectively.The phase φ(t) and amplitude A(t) are cal-culated individually at each point (x, y) respectively: where s(t) = H(s(t)) is the Hilbert transform of an arbitrary signal s(t), which is in our study either a mechanical signal derived from a deformation tensor or the time-course of the transmembrane potential.The phase angle φ(t) continuously increases over time from −π to π and uniquely represents the time-course of the action potential or a mechanical state, for instance, the deformation state within a cycle defined by two consecutive fully contracted states, see isoline for φ = −π in Fig. 2(c).As stated above, with isotropic tissues, we chose the mechanical signal to be the temporal evolution of the rotationally invariant quantity s(t) = Tr( )(t).In anisotropic tissues, we adjusted the mechanical signal to comprise only mechanical strains occurring along the fiber direction.For instance, with a uniform linearly transverse fiber orientation along the horizontal axis, we chose the signal to be s(t) = xx (t) axis correspondingly.In 3D anisotropic tissues, where the fiber direction rotates throughout the thickness of the bulk, we chose the signal to be the component of the Green-Lagrangian strain tensor corresponding to the local fiber axis, as later specified in section III.B. Prior to computing the phase φ and amplitude A, we subtracted the temporal average of the mechanical signal in each point: taken over the same time interval that the Hilbert transform is computed, where s is the original mechanical signal with a finite baseline.This procedure is used to determine the signal amplitude oscillating around its mean value at each point because time-averaged values of the signal can in general vary in space throughout the medium.Phase singularities can then be found in the phase maps φ(x, y) by computing the line integral [9]: where the integral is taken over a closed path and p and n are the numbers of positively and negatively charged phase singularities inside the area enclosed by the path, respectively.
In the anisotropic case, we additionally compute mechanical phase singularities by a different method used previously to identify the instantaneous tip position of spiral waves of electrical activity [32].In that context, the spiral wave shape was defined as a contour of constant transmembrane voltage U = U c , where U c is a constant chosen to lie in between the minimum and maximum values of U corresponding to the resting voltage and the         peak voltage following depolarization.This contour then separates resting and depolarized regions of the tissue and the spiral wave tip can be defined as the intersection of U = U c contour at time t and t + dt, or equivalently as the intersection of the U = U c and ∂ t U = 0 contours at time t.Similarly, we define the position of mechanical singularities as the intersection of the xx = c and ∂ t xx = 0 contours with the constant c defined as the average of the minimum and maximum values of xx .With this definition, the xx = c contour separates under-and over-strained regions of the tissue and mechanical singularities locate singular endpoints of this contour that can be interpreted as instantaneous deformation wave tips.The Hilbert transform and contour-intersection methods produce almost identical locations of mechanical singularities in most of the cases presented here where spatiotemporal patterns of voltage and mechanical activity are time-periodic.The contour intersection method is used here as additional validation of the Hilbert transform method as well as for physically interpreting the formation of unpaired mechanical phase singularities in terms of deformation wave breaks leading to the creation of deformation wave tips.
Figures 2(c-d) show the corresponding phase map φ and amplitude map A of the strain pattern shown in Fig. 2(b).The phase map in Fig. 2(c) reveals a pinwheel pattern, which retains a similar spiral-like shape, as in Fig. 2(a-b) and exhibits spiral-shaped lines of equal phase merging at the center of the medium.Using Eq. ( 17) one can detect a phase singular point (red dot) at the center of the medium, describing a topological defect point in the phase plane.At this point, the elastic medium may be contracted or dilated, but the signal does not change with time, and consequently, the amplitude approaches zero, see Fig. 2(d).Note that, due to the definition of amplitude (Eq.( 15)) in the Hilbert transformation, and the fact that we are working with the normalized response, amplitude being close to zero means that there is almost no difference between s and its 90 • shifted value (the difference between the real and imaginary parts of the Hilbert transformed signal).In other words, the signal remains almost constant in time.We will discuss this matter in more detail for the anisotropic spiral case.
What is clear from this simulation is that in the absence of fiber architecture (and consequently force anisotropy), the hydrostatic strain field Tr( ) follows the contraction field closely.As a result, the only mechanical phase singularity that we observe is due to the existing phase singularity of the prescribed active contraction force field.Therefore, we can hypothesize that in tissue with isotropic fiber architecture and no singularity of contraction field, we must not observe any mechanical phase singularity.To test this hypothesis, we apply an active contraction field of circular ring-shaped focal patterns to the tissue with isotropic fiber architecture.Again, the active contraction field was derived analytically as presented in Eq. ( 22) in Appendix A. The results of which are presented in Fig. 3.As one can see in Fig. 3, the contraction field T a (x, y) in (a), the strain field tr( )(x, y) or trace of the strain tensor in (b), and the phase map of the strain φ(Tr( ))(x, y) in (c) all retain the circular shape, and there is no phase singularity in the phase map.Moreover, the amplitude of the transformed signal remains relatively large everywhere inside the tissue.
In the cardiac muscle, the fiber architecture is highly anisotropic and organized in a complex, orthotropic pattern.Furthermore, the contracting active force is exerted only along the local fiber direction.As the next step, we study the effect of fiber anisotropy in the 2D tissue with a clockwise rotating spiral of active contraction.With no loss of generality, we assume that the fibers are aligned uniformly along the x-axis, and therefore, the contractile force occurs in this direction as well.The results for this case are presented in Fig. 4. Figure 4(a) shows the active contraction field T a (x, y).Note that the active contraction pattern is elongated along the horizontal x-axis, which reflects the underlying horizontal fiber anisotropy.The elongation occurs in real myocardium because the conduction is higher along the fiber direction,     and therefore, electrical current propagates faster in this direction.Since active contraction is slaved to electrical current, the active contraction will propagate faster and elongates along the same axis as well.Correspondingly, in the anisotropic simulation, the contractile force is only exerted along the horizontal x-axis, which represents myocardial muscle tissue more accurately than the isotropic simulations shown in Figs. 2 and 3. To identify mechanical phase singularities, while accounting for the effects caused by the fiber anisotropy, we analyzed the corresponding strain along the fiber orientation, which in this case is xx .Due to the effects of force anisotropy along the horizontal axis, the strain pattern in Fig. 4(b) is distorted compared to Fig. 2(b), and the generated strains are generally larger in the vertical arms.In other words, while the contraction field is similar in both vertical and horizontal spiral arms in both the aniso-and isotropic case, the strain is significantly diminished in the horizontal parts of the spiral arms in the anisotropic case because the tissue cannot generate sufficient contractile force perpendicularly to the fiber orientation.In comparison to the isotropic case, the phase map φ(x, y) in Fig. 4(c) is more complex, and it does not retain a simple pinwheel shape.The phase map reveals that, next to a mechanical phase singularity close to the core of the spiral pattern, additional mechanical phase singularities form at a distance from the central phase singularity in the horizontal spiral arms.The central mechanical phase singularity co-localizes with the rotational core or tip of the electrical (or active contraction) spiral pattern shown in Fig. 4(a), but the additional mechanical phase singularities do not co-localize with any phase singularities in the active contraction field (movies are available in SI).In this example, see Fig. 4(c), two pairs of additional mechanical phase singularities form in the vertical direction on each side of the core.The two singularities in each pair exhibit opposite topological charges, but do not exhibit vorticity as the central stationary phase singularity does.Overall, one notices an underlying pinwheel pattern, as in Fig. 2(b-c), but the pattern is topologically discontinuous in the space between each of the singularity pairs.Furthermore, the amplitude A in Fig. 4(d) vanishes in the vicinity of the mechanical phase singularities, independently of whether they are paired or unpaired with singularities of the contraction field.However, in the surrounding neighborhood of paired and unpaired mechanical phase singularities, the spatial distribution of the amplitude A differs significantly.Near paired singularities, the amplitude increases rapidly with distance from the singular point, while the amplitude vanishes at the center of the singular point in Fig. 4(d)).Near unpaired singularities, the amplitude remains small over a spatially extended region surrounding the singular point (black and red regions in Fig. 4(d)).As we show later, the property that the amplitude of oscillation grows rapidly or slowly away from paired or unpaired singularities, respectively, can be exploited to distinguish between those two different types of mechanical phase singularities.
To better understand the nature of unpaired mechanical phase singularities, we examine the strain's temporal evolution (measured along xx , here aligned with the fiber direction) close to and further away from an unpaired singularity, see Fig. 4(e).One can observe that in the vicinity of an unpaired mechanical singularity (e.g., point 1) the amplitude of xx vanishes over the entire period, whereas further away at a non-singular site xx exhibits large oscillations (e.g., point 2), c.f. Fig. 4(b).The inset in Fig. 4(e) shows how the corresponding trajectories of these two time-series look in a two-dimensional complex phase space with the real and imaginary parts of the signal on the x-and y-axis respectively: the non-vanishing time-series of the strain with finite amplitudes corresponds to a circular trajectory, whereas the strain with vanishing amplitude reduces to a point in phase space, which prohibits defining a phase angle φ( xx ).This behavior is rather unexpected for an unpaired singularity as shown in point 1 in Fig. 4(e).As we discuss next for the simpler case of a focal excitation, this property is the result of the combined effects of fiber anisotropy and the non-locality of elastic interactions that can cause xx to  remain constant in time.
The strain field resulting from a focal excitation wave in an anisotropic medium is shown in Fig. 5.The contraction field T a (x, y) and the strain field measured along the fiber axis xx (x, y) are shown in (a) and (b), respectively.Note that the thickness of the focal contraction wave is not uniform due to the differences in conduction speed of the excitation wave along the horizontal and vertical directions, respectively, similarly as in Fig. 4. As a result, the wave is thicker in regions in which it propagates in parallel to the fiber axis than in regions in which it propagates perpendicularly to that axis.As can be seen in Fig. 5(b), the strain in the vertically traveling parts of the elliptical wave is diminished.Importantly, panel (c) shows that the elliptical focal contraction wave produces four unpaired mechanical phase singularities in the tissue.These four phase singularities emerge in pairs and are aligned symmetrically around the focal center of the wave and travel along a confined closed loop (movies are available in SI).Lastly, panel (d) shows that the amplitude A of the transformed signal vanishes in the neigh-borhood of these phase singularities, as was the case with the anisotropic spiral wave pattern shown in Fig. 4.
To gain further insight into the mechanism of the formation of unpaired singularities, we paced an anisotropic tissue with a periodic wave train of electrical plane waves propagating from bottom to top along the y-direction while the fibers were aligned uniformly along the x-axis perpendicularly to the wave.A snapshot of this simulation is shown in Figs.6(a-b).While the active tension is spatially uniform parallel to the fiber axis in the depolarized contracting region of tissue, see Fig. 6(a), the resulting deformation wave measured by xx is "broken" in the middle region, see Fig. 6(b).This wave break gives rise to the formation of two mechanical singularities located at the tips of the xx = c contours (thin red dashed lines) separating under-strained and over-strained regions of tissue shown in green and red in Fig. 6(b), respectively.This wave break occurs because the active force in the central region of tissue is contracting against large regions of passive surrounding tissue, thereby producing small strains, while less constrained regions closer to the stress-free boundaries of the tissue produce a larger strain.In contrast, with fixed displacement boundary conditions on all edges of the tissue (u x = u y = 0), xx = 0 in the whole tissue is the trivial solution of linear elasticity under the same plane wave of active tension that does not produce any strain inhomogeneities.
However, unpaired mechanical singularities form independently of the particular mechanical boundary conditions with more complex wave patterns, such as the spiral and focal wave patterns.Figs.6(c-f liptical target wave propagating along the fiber axis are thicker and hence exert a larger net contractile force that overstrains those regions.In contrast, regions propagating perpendicular to this axis are thinner, thereby exerting a smaller net contractile force.Those regions are also more elongated, thereby contracting against a larger region of passive tissue.Both factors contribute to causing those thinner elongated regions to be understrained.A similar mechanism is seen to underlie the formation of unpaired mechanical singularities in the case of reentrant spiral waves, see Fig. 4, that, like elliptical target waves, exhibit both thicker and thinner regions propagating parallel and perpendicular to the fiber axis. We further investigated how the number of unpaired mechanical phase singularities increases with tissue size.Figure 7 shows the number of mechanical phase singularities per period as a function of ratio L/λ s of tissue size and spiral wavelength for a single spiral wave rotating in a 2D medium with anisotropic fiber architecture.The results show an increase in the number of mechanical phase singularities with system size.For small values of L/λ s (L/λ s = 1, 2) the only mechanical phase singularity appears at the core of the spiral, which is associated with the co-localized singularity at the tip of the spiral.As L/λ s = 3 the number of singularities jumps to 5 which accounts for the formation of two additional pairs of non co-localized singularities.This case (L/λ s = 3) is the same as the one presented in Fig. 4. As L/λ s increases further, the mechanical state of the medium becomes more complex and pairs of mechanical phase singularities are created and annihilated during each period, thereby causing the number of unpaired singularities to fluctuate during one period.However, the average number of phase singularities increases roughly linearly with the number of turns of the spiral wave that increases proportionally to L/λ s .This behavior is expected based on our physical interpretation of the formation of unpaired singularities based on deformation wave breaks.
In summary, our 2D simulations demonstrate that mechanical phase singularities co-localize with the tip of electrical spiral waves and characterize the center of an electromechanical rotor.In addition, unpaired mechanical phase singularities can form with both reentrant and focal waves, which then characterize the tip location of broken deformation waves.Wave breakup occurs due to the creation of overstrained and understrained regions of tissue that form preferentially along thicker and thinner regions of waves propagating parallel and perpendicular to the fiber axis, respectively.Since those waves can emanate from a focal source of excitation or a spiral wave center, the formation of unpaired mechanical singularity does not necessitate an electrical singularity to be present, thereby leading to a dissociation of a strict paired organization of electrical and mechanical phase singularities.In addition, mechanical boundary conditions can influence the formation of unpaired mechanical singularities, e.g., by facilitating their formation near stress-free surfaces during plane wave propagation, but unpaired singularities form independently of the boundary conditions in the case of reentrant or focal excitation waves.

B. Three-dimensional Dynamics
Phase singular points in 2D correspond to lines of phase singularity in 3D.Vortex dynamics in threedimensional excitable media have frequently been characterized using lines of phase singularity, where the lines or vortex filaments indicate the rotational core regions of three-dimensional scroll waves [32].Figure 8 depicts electromechanical phase singularities in a three-dimensional bulk tissue (1.2 × 1.2 × 1.2 cm 3 ) with (a) uniformly horizontally aligned fibers or linearly transverse fiber architecture and (b) fibers being rotated by 90 • along the z-axis, where at the bottom face of the bulk they are aligned along the x-axis and at the top surface they are aligned along the y-axis.The tissue is deformed by a single electrical scroll wave.The electrical scroll wave, which was simulated using the reaction-diffusion model presented in Sec.II B 1, can be described by a straight electrical vortex filament (black) being aligned vertically along the short axis in the center of the bulk in between its upper and lower surfaces.The electrical vortex filament does not bend or twist, and accordingly, the scroll wave remains stable.There are co-localized mechanical phase singularities (red filaments) for both cases shown The iso-surface is presenting Ēf = Ēxx = 0 where Ēf is the normalized strain along the fiber that is used to calculate the phase.
in (a) and (b).However, in the case of no fiber rotation, one extra "unpaired" mechanical phase singularity in the form of a filament attached to the boundary can be observed.In these simulations, the ratio between the size of the tissue and the spiral wavelength falls in the left part of the graph shown in Fig. 7. Therefore, there is only one pair of co-localized electrical and mechanical phase singularities describing the scroll wave core, if we ignore the one mechanical filament attached to the boundary.To study the effect of tissue size and wavelength on the spatial organization of mechanical filaments, we carried out a simulation with the same parameters as in Fig. 8, but in a larger tissue with size 4.5 × 4.5 × 1.2 cm, see Fig. 9.The filaments shown in Fig. 9 for a single scroll wave were obtained with no fiber rotation as in Fig. 8(a).A co-localized pair of an electrical and a mechanical filament mark the core region of the scroll wave at the center of the medium.Additionally, one can see that there are several unpaired mechanical phase singularities outside of the core region in analogy to the 2D anisotropic case shown in Fig. 4.
The excitation, contraction, and strain fields together with the phase, singularities, and amplitude on the top surface of the tissue are illustrated in Fig. 10.Panel (a) shows the electrical excitation U (dimensionless normalized units n.u.) on the top surface of the ventricular tissue slab intersected by the vortex filament.Due to the orientation of the scroll wave within the bulk, the near-surface electrical pattern corresponds to a clockwise rotating spiral wave pattern, and the tip of the spiral co- D K / w 5 q H 3 4 r 1 7 H 7 P W g j e P 8 B D + w P v 8 A f G K k X Y = < / l a t e x i t > FIG. 9. Electromechanical phase singularities in the presence of a single scroll wave with short wavelength in a 3D tissue (4.5 × 4.5 × 1.2cm, inset figure is the top view) with all fibers being aligned along the x-axis.Red lines represent lines of mechanical phase singularity or mechanical filaments and the black line at the center represents the electrical vortex filament or scroll wave core.The iso-surface is presenting Ēf = Ēxx = 0, where Ēf is the normalized strain along the fiber that is used to calculate the phase.
incides with the intersection point of the vortex filament with the surface.In addition, the electrical wave propagates faster along the x-axis parallel to the fibers (white arrow).This leads to a stretched spiral wave pattern in both electrical excitation U , as shown in (a), and active contraction T a , as shown in (b).To derive the strain field in (c) we computed the Green-Lagrangian strain tensor E as: where is the deformation gradient tensor, and I is the identity.Unlike in 2D, the 3D simulations exhibit large strains.The strain pattern was derived by computing the strain along the fiber, here denoted as E f , which in this case coincides with E xx , similarly as shown in Fig. 4(b).Note that the strain along the fiber i.e., E f , is a simple transformation of the strain tensor where the x−axis rotates to match the fiber direction.However, because the wavelength of the spiral is small in comparison to the medium size, the spiral is curled up more, and the anisotropic strain pattern can be observed not only in the near-field of the spiral wave core but for multiple wavelengths.The phase and amplitude maps     The spatial organization of unpaired singularities away from the central paired electromechanical singularity in Fig. 10 can be interpreted using the insights from Fig. 6 that highlighted the roles of 2D wave patterns (plane waves or focal excitations) and mechanical boundary conditions in the formation of unpaired singularities.The unpaired mechanical singularities forming far from the spiral core form preferentially near the stress-free surfaces of the tissue due to the breakup of nearly plane waves of deformation by a mechanism directly analogous to Fig. 6(a-b), while mechanical singularities closer to the spiral core region form by a mechanism analogous to the one illustrated in Fig. 6(c-d) for a focal excitation, which also pertains to reentrant waves.Furthermore, the migration of unpaired mechanical singularities closer to stress-free surfaces with distance from the spiral core in Fig. 10(d) can be interpreted as a transition between the near-core formation mechanism, which is insensitive to mechanical boundary conditions, and the far-core planewave breakup mechanism where stress-free surfaces favor the formation of unpaired singularities.
Finally, Fig. 11 depicts electromechanical vortex filament dynamics for the more realistic case that fibers are organized orthotropically and the fiber direction rotates along the z-axis through the thickness of the ventricular wall.In panels (a-f), with no loss of generality and for simpler visualizations, we show the normalized strain along the fiber Ēf .This measure was consequently used for further analysis using the Hilbert transform.At the bottom of the bulk at z = 0, the fibers are aligned uniformly along the x-axis.At the top, they are aligned along the y-axis.This results in a 90 • rotation over the thickness of the bulk of 1.2 cm.Panels (a-b) show the electromechanical filaments being produced by one clockwise rotating electrical scroll wave pattern.The comparison of Fig. 11(a-b) and Fig. 9 shows that the fiber rotation causes a more complex structure with bent filaments, but a qualitatively similar organization with a central paired electromechanical singularity (black and red) surrounded by additional unpaired mechanical singularities (red).
Next, panels (c-d) in Fig. 11 show the filament dynamics in the regime of electrical wave turbulence corresponding to steeper action potential duration restitution slope (Re = 1 instead of 0.8, see Eq. ( 30)).In this regime, scroll wave breakup leads to a multiplication of electrical filaments.Interestingly, the results in Figs.11(c-d) show that the electrical (black) and mechanical filaments (red) are not as well co-localized in the scroll wave breakup regime compared to the single stable scroll wave regime of Figs.11(a-b).However, paired electromechanical singularities can still be distinguished from unpaired singularities that are not co-localized with any electrical filaments.Lastly, we show in Figs.11(e-f) the results obtained with a 3D focal wave after pacing the bulk at the center of its top surface.The results confirm that unpaired mechanical singularities can exist without the presence of electrical singularities as a direct extension of the 2D results shown in Fig 5 .To conclude the results section, we present a method to distinguish between paired electromechanical singularities and unpaired mechanical singularities, which we base on our observations of the behavior of the complex amplitude maps surrounding the mechanical singularities, see Fig. 4(d) and Fig. 10(e).We exploit that the amplitude of deformation waves is significantly smaller in the vicinity of unpaired mechanical singularities compared to paired electromechanical phase singularities.To do so in a consistent way, we first normalize the amplitude map between 0 and 1 at each time step.Then ¸A(r)dl is calculated around a mechanical phase singular point, where A is the amplitude and l is a closed path around that point.This closed path is one grid square, i.e., ˛Ai,j (r)dl = A i,j + A i+1,j + A i+1,j+1 A singular point is then selected as being paired to an electrical singularity or unpaired depending on whether the amplitude is larger or smaller than some threshold In the case that electrical waves remain stable one can see a clear co-localization between mechanical and electrical filaments.However, as scroll waves break, a more complex mechanism can be observed, but many of electrical filaments are still paired with a mechanical one.
(e.g., 0.15).Fig. 12 shows how paired and unpaired mechanical phase singularities can be discriminated using this filtering method.The figure shows the spatiotemporal organization of electrical and mechanical phase singularities on the top surface of the bulk tissue.Panel (a) depicts the case of a single spiral and linear transverse fiber alignment or no fiber rotation, c.f. Fig. 9. Panel (b) presents the same case, but with fiber rotation, c.f. Fig. 11(a).Lastly, panel (c) depicts the case of multiple scroll waves due to electrical wave break in the tissue with fiber rotation, c.f. Fig. 11(b).In this figure, the electrical and mechanical phase singularities emerge on the surface of the bulk tissue (sampled over a short time interval ≈ 0.4s), and are shown as black circles and red crosses, respectively, and the filtered mechanical phase singularities are shown as blue circles.One can see that in cases (a) and (b), where there is a single scroll wave and no electrical wave break, the filtering method accurately selects the paired mechanical phase singularities.In the wave turbulence regime, the method is still reasonably effective at eliminating the unpaired mechanical phase singularities and the reduced accuracy can be attributed to the aforementioned weaker co-localization of paired electrical and mechanical singularities in this regime.In summary, it is possible to automatically distinguish paired from unpaired mechanical singularities using only information about the temporal evolution of the strain in proximity to a singularity.

IV. DISCUSSION
In this study, we showed that electromechanical phase singularities are an integral and universal phenomenon in elastic excitable media and result from the dynamics and coupling between cardiac excitation waves, anisotropic contraction, and tissue strain.Electromechanical phase singularities are composed of electrical and mechanical phase singularities, which emerge in the tissue's electrophysiology and elasticity, respectively.They describe similar features of the respective dynamics and have a distinct relationship to each other while exhibiting separate characteristics pertaining to the particular physics of excitation waves and tissue strain, respectively.While we found that electrical and mechanical phase singularities can emerge as co-localized pairs, in which case they equally describe the core of an electrical spiral or scroll wave, we also found that mechanical phase singularities can furthermore form independently from electrical ones.
Our study provides for the first time a systematic and detailed description of the behavior of electromechanical phase singularities in a range of situations, in 2D and 3D media with and without muscle fiber anisotropy and different boundary conditions, with simple and more complex dynamics, such as focal waves, spiral waves and scroll wave chaos, and with fully relaxed analytically derived strain fields versus strain fields obtained in reactiondiffusion-mechanics simulations.We were particularly interested in the properties of mechanical phase singularities and in the questions of whether they (i) follow a similar organization as electrical ones and (ii) how re-liably they can be used to locate the core region of the electrical spiral and scroll waves that underlie cardiac fibrillation [9,10,14].
Our study demonstrates that mechanical phase singularities are topological defects, which emerge in the strain dynamics of contracting cardiac tissue due to both rotational electrical excitation waves and muscle fiber anisotropy.Because active contraction is exerted along muscle fibers, the deformation resulting from a particular excitation wave pattern is highly anisotropic, which manifests as a polarized strain field when measuring tissue strain.The anisotropy and polarization lead to mechanical wavebreaks, which are characterized by mechanical singularities.Importantly, mechanical singularities can also emerge in electrical regimes, which do not exhibit electrical phase singularities: focal excitation waves can produce mechanical phase singularities in the presence of anisotropy.We should emphasize that in this article we focused on elucidating the mechanism for the creation of mechanical singularities using minimal representations of the electrophysiology and mechanical response.This allowed us to rigorously dissect the problem.There is no question that our approach neglects the complex heterogeneous structure and geometry of the anatomical heart.However, we believe that our extensive set of simulations clearly shows that mechanical wavebreaks are the origins of these singularities.For instance, in the linearly transverse anisotropic case, an elliptical ring-shaped excitation wave produces two mechanical waves propagating in opposite directions along the fiber direction away from the focus, while the strain vanishes in the perpendicular direction to the fiber direction, see Figs. 5(b)  and 6(d,f).The elliptical excitation wave and the two resulting mechanical waves are characterized by four mechanical phase singularities, organized in two pairs surrounding the focus, where each pair describes one of the mechanical waves, and the phase singularities indicate the mechanical wavebreak locations, respectively.
In the generic situations that we studied in Figs.2-5, electromechanical phase singularities can be divided into paired or unpaired electrical and mechanical singularities.However, in more complex situations, their relationship can be more complicated.Paired electrical and mechanical singularities are co-localized, only exist in the presence of spiral and scroll waves, and emerge equally with and without anisotropy, c.f. Figs. 2 and 4. Paired mechanical singularities originate from the rotation of a strain field that is enslaved to the rotation of an electrical excitation wave, where the pairing occurs close to the tip or core region of the electrical spiral or scroll wave, which is typically characterized by an electrical singularity.Unpaired mechanical and electrical singularities occur in the presence of anisotropy with both focal and spiral or scroll waves because, in anisotropic tissues, focal waves inherently produce mechanical singularities, and spiral or scroll waves produce additional unpaired mechanical singularities outside of their core regions, as shown in Figs.4(b,c), 8 and 9.However, during three-dimensional chaotic scroll wave dynamics, we observed neither purely paired nor unpaired, but instead partially paired electromechanical filaments close to the scroll wave's core.The degree of co-alignment of the partially paired filaments varies in space and over time, respectively; see Fig. 11(c,d).At this point, we can only speculate about the origins of this phenomenon, it could be that scroll wave drift or the complex fiber rotation cause the partial dissociation, or it could be that the local strain field is distorted by neighboring waves and their contractions, and the phenomenon needs further investigation.
From an electrophysiological perspective, mechanical phase singularities may carelessly be disregarded as an ambiguous, even unreliable concept to describe electrical rotors because unpaired mechanical singularities do not co-localize with electrical singularities and could accordingly be interpreted as false positives.Further, paired mechanical singularities are merely proxies for the positions of electrical singularities.However, while we demonstrated that it is possible to automatically distinguish paired from unpaired mechanical singularities, see Fig. 12, therefore demonstrating that the localization of electrical singularities via paired mechanical singularities would be practically feasible, our research also shows that mechanical singularities are much more than a mere description of electrical phenomena.They constitute a more generalized way to characterize electromechanical tissue dynamics because they simultaneously reflect the topology of excitation wave patterns, strain dynamics, and interactions of the excitation with the underlying mechanical substrate.The spatiotemporal organization of tissue deformation is far more complex than the spatiotemporal organization of electrical waves, and mechanical phase singularities reflect this more complex organization.This is exemplified by spiral and focal excitation waves, each producing a very characteristic strain and mechanical singularity pattern.Mechanical singularities encode both the principal morphology of the excitation wave pattern as well as the underlying organization of muscle fibers and the deformation state of the muscle.For instance, single high-frequency focal sources produce a complex pattern of unpaired mechanical filaments, see Fig. 11(e-f), and the filaments are a signature of this particular activation pattern together with a specific underlying muscle fiber organization, and not just spurious topological defects.To further avoid the impression that mechanical singularities are poor proxies for electrical singularities, we would like to emphasize that unpaired mechanical singularities are only numerous in a regime exemplified in Figs. 9 and 10 where the tissue size far exceeds the spiral wavelength, which is uncommon during ventricular fibrillation.Figs. 4 and 8 are more representative of an episode of polymorphic ventricular tachycardia or fibrillation with fewer unpaired mechanical singularities when the wavelength of reentrant waves is comparable to the tissue size.Therefore, we expect the number of unpaired singularities to be typically comparable to the number of paired singularities during high-frequency cardiac arrhythmias.Unpaired singularities may even be completely absent for slower forms of anatomical tachycardias with large wavelengths.
One limitation of the present study is that it rests on the assumption that the calcium transient, and hence mechanical contraction, is in phase with electrical excitation.This assumption remains valid as long as calcium release from intracellular calcium stores is triggered by calcium entry into the cell via L-type calcium channels following membrane depolarization.The fact that spiral waves can be imaged in heart tissue from an optical mapping of the calcium signal [36] provides direct evidence that this calcium-induced-calcium-release (CICR) mechanism, which underlies normal physiological function [15], can remain operative during reentrant cardiac arrhythmias.However, in other settings, the calcium transient may be suppressed if the activation interval is too short for intracellular stores to refill with calcium during the action potential, or may no longer be in phase with the voltage signal, e.g. if the calcium transient occurs by spontaneous release from intracellular stores instead of CICR following membrane depolarization.In those settings, even paired mechanical singularities may cease to exist.Another limitation of the present study is that it neglects mechanoelectrical feedback [16].This feedback, which is complex and still not completely understood, originates from the fact that stretch can modify passive constitutive properties such as membrane capacitance and electrical coupling and activate mechanosensitive ion channels.This feedback has been shown to influence properties of excitation waves such as action potential duration restitution, conduction velocity, and to even induce spiral drift [18].However, even in the presence of mechanoelectrical feedback, the wave pattern of active force still tracks closely the one of electrical excitation as long as CICR remains operative, e.g. a drifting spiral wave of excitation will induce a drifting pattern of active force.The resulting field of tissue deformation, however, will still exhibit a different spatiotemporal organization including wavebreaks due to the nonlocality of elastic interactions.Therefore, we do not expect that mechanoelectrical feedback would qualitatively change our main finding of the existence of both paired and unpaired mechanical singularities.
In the future, our findings could play an important role in the interpretation of ultrasound imaging data of the fibrillating heart.Mechanical phase singularities could be used to characterize the three-dimensional tissue dynamics during ventricular or atrial fibrillation and could provide estimates of the anchoring sites of electrical spiral or scroll waves deep within the tissue.To further assess the feasibility of a mechanics-based arrhythmia imaging technique, we aim to study similar electromechanical phenomena beyond the generic setting in physiologically detailed models, explore potential limitations and determine the effect of heterogeneity and tissue anatomy, as mechanics will be substantially altered in fibrotic or scar tissue.Lastly, we want to better understand the dissociated electromechanical filament dynamics encountered during scroll wave chaos, which is likely caused by the long-range character of elastic forces.This requires addressing the fundamental question of whether the inverse mapping between the observed deformation wave pattern and the excitation-contraction wave pattern causing the deformation, is unique.In recent work it was shown in silico that it is possible to compute even complicated electrical excitation wave patterns, such as chaotic scroll waves, from tissue deformation [34,37], suggesting that it could also be possible to find a unique mapping between electrical and mechanical filaments.

ACKNOWLEDGMENTS
The authors acknowledge support from the Center for Interdisciplinary Research on Complex Systems at Northeastern University.J.C. and S.L. acknowledge support from the German Center for Cardiovascular Research, partnersite Göttingen.J.C. acknowledges support from the Cardiovascular Research Institute at the University of California, San Francisco.S.L. acknowledges support from the Max Planck Society.The simulations benefited from computing time allocation on Northeastern University Discovery Cluster at the Massachusetts Green High Performance Computing Center (MGHPCC).

A. ANALYTICAL EXPRESSIONS FOR SPIRAL AND TARGET WAVES
The active contraction field, T a for an n-turn Archimedean spiral, rotating at angular frequency ω for isotropic conduction is given by where θ(x, y, t) = arctan −x sin ωt + y cos ωt x cos ωt + y sin ωt .
λ s = (2πc)/ω is the spiral wavelength and c is the conduction velocity.σ s represents the thickness of the spiral arms.
For the pacing case, the active contraction field is shaped as a target wave pattern generated by a sequence of n stimuli delivered at (x = 0, y = 0) with period T .The contraction field of the isotropic pacing case is: where σ s is the wave thickness as in the spiral wave case, r i (t) = c(t − iT )Θ(t − iT ) is the radius of the i th concentric wave, c is the conduction velocity, and Θ(x) is the Heaviside step function defined by Θ(x) = 1 for x ≥ 0 and Θ(x) = 0 for x < 0.
The anisotropic case where fibers are aligned along the x−axis can simply be achieved by rescaling the x coordinate (x → x/ D /D ⊥ ) in Eqs.(20) and (22). in this rescaling D /D ⊥ is the ratio of the diffusivity parallel and perpendicular to the fiber direction, chosen equal to √ 5 in our simulations.

B. REACTION-DIFFUSION MODEL
The main reaction-diffusion equations to model the electric field has the form: where Re is the parameter that controls the restitution properties.Increasing Re makes the slope of the action potential duration restitution curve steeper at the short diastolic interval, thereby promoting spiral/scroll wave breakup.The parameter M controls the conduction velocity restitution curve and increasing M flattens this curve.In all of the results that are presented in this article except one case, Re = 0.8 and M = 10, which corresponds to a regime where a single spiral or scroll wave propagates and does not break up.The only exception is the investigation of a 3D parallelepipedal slab of tissue with rotation anisotropy where Re = 1.0 produces scroll wave breakups.More details about this model and the effect of each parameter can be found in [26].In Eq. ( 7), χ(U ) is a step function that sets the time scale of the contraction period with respect to U .It has been set as χ(U ) = 5 if U ≥ U h otherwise, χ(U ) = 20.The set of parameters that are used in this article are listed in Table I.A characteristic pulse and contraction structure based on these parameters are shown in Fig. 13.In this figure, a 1D discretization of Eqs.( 5), (6), and (7) has been solved using the forward Euler method.One can see that during a single pulse, the voltage surges immediately to its maximum.However, the contracting force starts to develop slower than the voltage, and it reaches its maximum when the voltage almost starts to decrease.Even though this model is phenomenological, it reproduces qualitatively the delayed peak contractile force following depolarization.

C. ANALYSIS OF 1D CABLE
In a human heart with density ρ ≈ 1000 kgm −3 and Young's modulus E ≈ 125 MPa [18], the elastic wave speed can be crudely estimated to be C ≈ 353 ms −1 .Furthermore, the period of excitation T spans the range of 0.1 to 1.0 seconds with the lower and upper bounds corresponding to an electrical rotor and a normal heartbeat, respectively.It follows that the distance CT traveled by elastic waves during a contraction period is much larger than the characteristic size L of the heart of a few centimeters.This scale separation (L CT ) makes it possible to rescale the equations for the computational efficiency of the time integration algorithm.I.The left vertical axis corresponds to the variables U or v and the right vertical axis represents the dimensionless magnitude of the active contracting strain Ta.The horizontal axis is time in the units of milliseconds.

A. Spatial decay rate of elastic waves
To choose our model parameters to be in a regime where L CT , we compute the spatial decay rate of elastic waves starting from the 1D elastodynamics equation with Kelvin-Voigt damping and no external force or equivalently Before computing the decay rate, it is useful to first derive expressions for the storage and loss moduli in a setting where a one-dimensional cable is held fixed at one end (x = 0) and periodically displaced at the other end (x = L) corresponding to the boundary condition u(0) = 0; u(L) = A cos(ωt) (34) where E and η are Young's modulus and damping of the system, C = E/ρ is the characteristic wave speed of the system, D e = η/ρ is the elastic diffusion constant of the system and L is the length of the cable.With L CT (T = 2π/ω), Eq. (32) will reduced to E∂ 2 x u + η∂ 2 x u = 0 (35) One solution that satisfies this equation has the form of which, together with ρ∂ 2 t u = ∂ x σ, yields σ = E∂ x u + η∂ x u = EA L {cos(ωt) − τ R ω sin(ωt)} (37) where τ R = η/E is the retardation time.Further using the trigonometry identity cos(ωt + δ) = cos(ωt) − δ sin(ωt) in the limit of small δ we will have where δ = ηω/E and σ 0 = EA/L.Using Eq. ( 36), the strain in the cable is where 0 = A/L.In Eq. (38), δ is the phase lag between stress and strain.Using Eqs. ( 38) and (39) the storage and loss moduli are Next, to calculate the spatial decay rate of elastic waves, we substitute a traveling wave solution of the form into the full elastodynamics equation including inertia (Eq.( 32)), which yields or where the magnitude of the imaginary part of k determines the decay rate B. Choice of model parameters In our simulations, we set C = 5 ms −1 and the size of the tissue is 4.5 × 4.5 × 1.2 cm 3 .Therefore, in the limits of 1 s and 0.1 s periods, we have: Those estimates are consistent with our assumption that that the distance CT traveled by elastic waves during a contraction period is much larger than the characteristic size L of the heart of a few centimeters.

C. Adimensionalization
In this section, we adimensionalize Eq. ( 33) and Eq. ( 13) of the MSM using a characteristic length L 0 = a where a is the lattice size, a characteristic time, τ el = L 0 /C, and ρ.Substitute these characteristic values into Eq.( 33) and simplifying, one will have: In this fashion, the dimensionless characteristic wave speed of the system is unity, and the dimensionless damping of the system is D e /(Ca).The corresponding discrete elastodynamics equation for the 1D MSM is where F is the sum of active and volumetric forces applied to the mass.Let us adimensionalize Eq. (49) (the same adimensionalization holds in 3D for Eq. ( 13)) by the characteristic length of the lattice size a, E the Young's modulus of the tissue, τ el the characteristic time of the system, and ρ the density of the material, yielding To characterize the mechanical properties of the massspring model, such as stress-strain relation or storage and loss modulus, we conducted a series of numerical studies.To approximate one-dimensional settings, we model a long bar on a cubic lattice where all the simulation parameters are presented in Table II.We used ρ, E, and a to adimensionalize the equation of motion utilizing equations presented in part C of Appendix C.  As shown in Fig. 14(a), a pair of tensile displacements along x−axis has been applied on both ends of the bar with traction free boundary conditions on the top and the bottom surfaces.To calculate Young's modulus and Poisson ratio, for a given imposed strain, we obtain the quasistatic configuration by letting velocities and accelerations of all masses to vanish.To speed the simulations, we set the adimensionalized damping constant D e /(Ca) = 100 where D e = η/ρ is the elastic diffusivity and η is the damping ratio.Note that the small time step is because of the large damping we used for the simulations.
First, we investigated the effect of the volume conservation penalty pressure on the elastic properties of the MSM.The analytic equations for these relations for infinitesimal strains have been presented in Eqs.11 and 12.We studied these properties of the MSM for strains xx = u x /L = 0.001, 0.15.The former represents the infinitesimal strain elasticity limit and the latter can be considered as the large deformation limit.The 0.15 strain limit was chosen because, during a traveling electrical wave, cardiac cells experience about 15% compression strain along their fibers [18].In these simulations, stresses are calculated by adding all the reaction forces at one end and dividing it by the engineering crosssection (5a × 5a).Using the calculated stress and applied strain, the Young's Modulus of the MSM is calculated as E = σ xx / xx .In addition, we calculated the Poisson's ratio buy finding the ratio between the transverse strain at the middle of the bar and the applied longitudinal strain i.e., ν = − yy / xx .Results as a function of the volumetric pressure are presented in Fig. 15.
In the small strain limit = 0.1%, the simulation results follow the analytical equations of 11 and 12 perfectly.However, for large strain cases, due to geometric nonlinearities, the MSM exhibits a hyperelastic behavior with higher Young's modulus and lower Poisson's ratio.In our simulations for the heart tissue, we assume p = p/E = 3 that yields the Poisson's ratio range be-tween 0.41 and 0.45 for the strain range between 0.1% and 15%.
To study the stress-strain relation of the MSM, we incrementally increased the applied displacement (strain) while volumetric penalty pressure was kept constant at a given values 0 ≤ p ≤ 3. The results, depicted in Fig. 16, clearly show that the presented MSM has a hyperelastic response, where increasing the volumetric penalty pressure will increase the nonlinearity of the model.Finally, in the small deformation limit we studied the viscoelastic properties of the MSM.We applied a sinusoidal displacement (strain) with 1 Hz frequency at x = L while the degrees of freedom at x = 0 were fixed.We then calculated the stresses from the reaction forces at x = L (see Fig. 14b).We define the mechanical transfer function (the complex modulus of the system) as Ẽ = σxx /˜ xx where σxx and ˜ xx are Fourier transformed axial stress and strain, respectively.We then define the storage modulus as the real part E s = Re( Ẽ) and the loss modulus as the imaginary part E l = Im( Ẽ) of the complex modulus.In a viscoelastic material, the ratio between the loss and storage modulus is a measure of damping in the material.This ratio as a function of dimensionless elastic diffusion is shown in Fig 17 .The elastic diffusion can be defined as the damping coefficient divided by the density of the material (see Appendix C part C).As expected, with zero damping there is no loss modulus, and the ratio between loss and storage modulus change linearly with the damping of the system.The effect of volumetric penalty pressure is not significant.Details are discussed in Appendix C. FIG. 17.The ratio between loss and storage modulus as a function of dimension less elastic diffusion.De is η/ρ where η is the damping of the system and ρ is the density of the material.C = E/ρ is the characteristic elastic wave speed and L is the length of the system (not the lattice dimension).The top x−axis shows the ratio between De and the diffusion coefficient D = 1.1 cm 2 /s in the electrophysiology model.The relation between E l /Es and damping of the system is linear in 0.1% strain and at zero damping, loss modulus is zero.The penalty volumetric pressure p slightly decrease the E l /Es ratio in high damping constants.
t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " > A A A C E H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u C G 5 c V n b b Q D i W T Z t r Q J D M k m U I Z + g l u X O i v u B O 3 / o F / 4 t J M O w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n 6 2 r u m f 5 w 3 W t 0 S l K K 8 M Z n M M l e H A D D b i H J v h A Y A j P 8 A p v z o v z 7 n w 4 n 4 v R k l P s n M I S n K 9 f A G q d 8 g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " > A A A C E H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u C G 5 c V n b b Q D i W T Z t r Q J D M k m U I Z + g l u X O i v u B O 3 / o F / 4 t J M O w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n 6 2 r u m f 5 w 3 W t 0 S l K K 8 M Z n M M l e H A D D b i H J v h A Y A j P 8 A p v z o v z 7 n w 4 n 4 v R k l P s n M I S n K 9 f A G q d 8 g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " > A A A C E H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u C G 5 c V n b b Q D i W T Z t r Q J D M k m U I Z + g l u X O i v u B O 3 / o F / 4 t J M O w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n 6 2 r u m f 5 w 3 W t 0 S l K K 8 M Z n M M l e H A D D b i H J v h A Y A j P 8 A p v z o v z 7 n w 4 n 4 v R k l P s n M I S n K 9 f A G q d 8 g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " > A A A C E H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u C G 5 c V n b b Q D i W T Z t r Q J D M k m U I Z + g l u X O i v u B O 3 / o F / 4 t J M O w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n

4 <
6 2 r u m f 5 w 3 W t 0 S l K K 8 M Z n M M l e H A D D b i H J v h A Y A j P 8 A p v z o v z 7 n w 4 n 4 v R k l P s n M I S n K 9 f A G q d 8 g = = < / l a t e x i t > c l a t e x i t s h a 1 _ b a s e 6 4 = " C d g q S N 7 P 6 8 L l Z G d 0 r Q 5 z 6 D G 0 B F 8 = " > A A A C G H i c b V C 7 S g N B F L 0 b X z G + o p Y 2 g 0 G w C r s i a B m w s Y x g H p A s Y X Y y m w y Z x z I z K 4 Q l n 2 F j o b 9 i J 7 Z 2 / o m l s 8 k W J v H A w O G c e 7 l n T p R w Z q z v f 3 u l j c 2 t 7 Z 3 y b r m i T W y W 7 m e R W M l j 4 9 s w Y z J J L Z V k E S d O O b I K 5 S 2 h I d O U W D 5 1 B B P N 3 I 8 Q G W P X l X V d V l x V w W o x 6 6 R 9 V Q 8 c f 7 i u N b p F a W U 4 g 3 O 4 h A B u o A H 3 0 I Q W E F D w D K / w 5 r 1 4 7 9 6 H 9 7 k Y L X n F z i k s w f v 6 B a 2 0 o Y Y = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " C d g q S N 7 P 6 8 L l Z G d 0 r Q 5 z 6 D G 0 B F 8 = " > A A A C G H i c b V C 7 S g N B F L 0 b X z G + o p Y 2 g 0 G w C r s i a B m w s Y x g H p A s Y X Y y m w y Z x z I z K 4 Q l n 2 F j o b 9 i J 7 Z 2 / o m l s 8 k W J v H A w O G c e 7 l n T p R w Z q z v f 3 u l j c 2 t 7 Z 3 y b r m i T W y W 7 m e R W M l j 4 9 s w Y z J J L Z V k E S d O O b I K 5 S 2 h I d O U W D 5 1 B B P N 3 I 8 Q G W P X l X V d V l x V w W o x 6 6 R 9 V Q 8 c f 7 i u N b p F a W U 4 g 3 O 4 h A B u o A H 3 0 I Q W E F D w D K / w 5 r 1 4 7 9 6 H 9 7 k Y L X n F z i k s w f v 6 B a 2 0 o Y Y = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " C d g q S N 7 P 6 8 L l Z G d 0 r Q 5 z 6 D G 0 B F 8 = " > A A A C G H i c b V C 7 S g N B F L 0 b X z G + o p Y 2 g 0 G w C r s i a B m w s Y x g H p A s Y X Y y m w y Z x z I z K 4 Q l n 2 F j o b 9 i J 7 Z 2 / o m l s 8 k W J v H A w O G c e 7 l n T p R w Z q z v f 3 u l j c 2 t 7 Z 3 y b r m i T W y W 7 m e R W M l j 4 9 s w Y z J J L Z V k E S d O O b I K 5 S 2 h I d O U W D 5 1 B B P N 3 I 8 Q G W P X l X V d V l x V w W o x 6 6 R 9 V Q 8 c f 7 i u N b p F a W U 4 g 3 O 4 h A B u o A H 3 0 I Q W E F D w D K / w 5 r 1 4 7 9 6 H 9 7 k Y L X n F z i k s w f v 6 B a 2 0 o Y Y = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " C d g q S N 7 P 6 8 L l Z G d 0 r Q 5 z 6 D G 0 B F 8 = " > A A A C G H i c b V C 7 S g N B F L 0 b X z G + o p Y 2 g 0 G w C r s i a B m w s Y x g H p A s Y X Y y m w y Z x z I z K 4 Q l n 2 F j o b 9 i J 7 Z 2 / o m l s 8 k W J v H A w O G c e 7 l n T p R w Z q z v f 3 u l j c 2 t 7 Z 3 y b

2 <
r m i T W y W 7 m e R W M l j 4 9 s w Y z J J L Z V k E S d O O b I K 5 S 2 h I d O U W D 5 1 B B P N 3 I 8 Q G W P X l X V d V l x V w W o x 6 6 R 9 V Q 8 c f 7 i u N b p F a W U 4 g 3 O 4 h A B u o A H 3 0 I Q W E F D w D K / w 5 r 1 4 7 9 6 H 9 7 k Y L X n F z i k s w f v 6 B a 2 0 o Y Y = < / l a t e x i t > c l a t e x i t s h a 1 _ b a s e 6 4 = " H e F S F c / L L U 6 U W e N i o m i B U A / i T y I X e v B f v 3 f v w P h e j G 1 6 x c w Z L 8 L 5 + A a p e o Y Q = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " H e F S F c / L L U 6 U W e N i o m i B U A / i T y I X e v B f v 3 f v w P h e j G 1 6 x c w Z L 8 L 5 + A a p e o Y Q = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " H e F S F c / L L U 6 U W e N i o m i B U A / i T y I X e v B f v 3 f v w P h e j G 1 6 x c w Z L 8 L 5 + A a p e o Y Q = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " H e F S F c / L L U 6 U W e N i o m i B U A / i T y I

2 <
X e v B f v 3 f v w P h e j G 1 6 x c w Z L 8 L 5 + A a p e o Y Q = < / l a t e x i t > k l a t e x i t s h a 1 _ b a s e 6 4 = " p e b E T C x I l V M F 2 P C m V y h P T n d s r m i T W y W 7 m e R W M l j 4 9 s w Y z J J L Z V k E S d O O b I K 5 S 2 h I d O U W D 5 1 B B P N 3 I 8 Q G W P X l X V d l l 1 V w W o x 6 6 R d r w W O P 1 x X G 9 2 i t B K c w w V c Q Q A 3 0 I B 7 a E I L C C h 4 h l d 4 8 1 6 8 d + / D + 1 y M b n j F z h k s w f v 6 B b f O o Y w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " p e b E T C x I l V M F 2 P C m V y h P T n d s r m i T W y W 7 m e R W M l j 4 9 s w Y z J J L Z V k E S d O O b I K 5 S 2 h I d O U W D 5 1 B B P N 3 I 8 Q G W P X l X V d l l 1 V w W o x 6 6 R d r w W O P 1 x X G 9 2 i t B K c w w V c Q Q A 3 0 I B 7 a E I L C C h 4 h l d 4 8 1 6 8 d + / D + 1 y M b n j F z h k s w f v 6 B b f O o Y w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " p e b E T C x I l V M F 2 P C m V y h P T n d s r m i T W y W 7 m e R W M l j 4 9 s w Y z J J L Z V k E S d O O b I K 5 S 2 h I d O U W D 5 1 B B P N 3 I 8 Q G W P X l X V d l l 1 V w W o x 6 6 R d r w W O P 1 x X G 9 2 i t B K c w w V c Q Q A 3 0 I B 7 a E I L C C h 4 h l d 4 8 1 6 8 d + / D + 1 y M b n j F z h k s w f v 6 B b f O o Y w = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " p e b E T C x I l V M F 2 P C m V y h P T n d s

4 <
r m i T W y W 7 m e R W M l j 4 9 s w Y z J J L Z V k E S d O O b I K 5 S 2 h I d O U W D 5 1 B B P N 3 I 8 Q G W P X l X V d l l 1 V w W o x 6 6 R d r w W O P 1 x X G 9 2 i t B K c w w V c Q Q A 3 0 I B 7 a E I L C C h 4 h l d 4 8 1 6 8 d + / D + 1 y M b n j F z h k s w f v 6 B b f O o Y w = < / l a t e x i t > k l a t e x i t s h a 1 _ b a s e 6 4 = " W B + k X 6 d R n 1 E B 4 r i j q B 8 c Y z c E + u I = " > A A A C G H i c b V C 7 S g N B F L 0 b X z G + o p Y 2 g 0 G w C r s i a B m w s Y x g H p A s Y X Y y m w y Z x z I z K 4 Q l n 2 F j o b 9 i J 7 Z 2 / o m l s 8 k W J v H A w O G c e 7 l n T p R w Z q z v f 3 u l j c 2 t 7 Z 3 y b m V v / + D w q H p 8 0 j Y q 1 Y S 2 i O J K d y N s K G e S t i y z n H Y T T b G I O O 1 E k 7 v c 7 z x R b Z i S j 3 a a 0 F D g k W Q x I 9 g 6 q d e P N S b Z Z J Z d z w b V m l / 3 5 0 D r J C h I D 6 e w B O / r F 7 s k o Y 4 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " W B + k X 6 d R n 1 E B 4 r i j q B 8 c Y z c E + u I = " > A A A C G H i c b V C 7 S g N B F L 0 b X z G + o p Y 2 g 0 G w C r s i a B m w s Y x g H p A s Y X Y y m w y Z x z I z K 4 Q l n 2 F j o b 9 i J 7 Z 2 / o m l s 8 k W J v H A w O G c e 7 l n T p R w Z q z v f 3 u l j c 2 t 7 Z 3 y b m V v / + D w q H p 8 0 j Y q 1 Y S 2 i O J K d y N s K G e S t i y z n H Y T T b G I O O 1 E k 7 v c 7 z x R b Z i S j 3 a a 0 F D g k W Q x I 9 g 6 q d e P N S b Z Z J Z d z w b V m l / 3 5 0 D r J C h I D 6 e w B O / r F 7 s k o Y 4 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " W B + k X 6 d R n 1 E B 4 r i j q B 8 c Y z c E + u I = " > A A A C G H i c b V C 7 S g N B F L 0 b X z G + o p Y 2 g 0 G w C r s i a B m w s Y x g H p A s Y X Y y m w y Z x z I z K 4 Q l n 2 F j o b 9 i J 7 Z 2 / o m l s 8 k W J v H A w O G c e 7 l n T p R w Z q z v f 3 u l j c 2 t 7 Z 3 y b m V v / + D w q H p 8 0 j Y q 1 Y S 2 i O J K d y N s K G e S t i y z n H Y T T b G I O O 1 E k 7 v c 7 z x R b Z i S j 3 a a 0 F D g k W Q x I 9 g 6 q d e P N S b Z Z J Z d z w b V m l / 3 5 0 D r J C h I D 6 e w B O / r F 7 s k o Y 4 = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " W B + k X 6 d R n 1 E B 4 r i j q B 8 c Y z c E + u I = " > A A A C G H i c b V C 7 S g N B F L 0 b X z G + o p Y 2 g 0 G w C r s i a B m w s Y x g H p A s Y X Y y m w y Z x z I z K 4 Q l n 2 F j o b 9 i J 7 Z 2 / o m l s 8 k W J v H A w O G c e 7 l n T p R w Z q z v f 3 u l j c 2 t 7 Z 3 y b m V v / + D w q H p 8 0 j Y q 1 Y S 2 i O J K d y N s K G e S t i y z n H Y T T b G I O O 1 E k 7 v c 7 z x R b Z i S j 3 a a 0 F D g k W Q x I 9 g 6 q d e P N S b Z Z J Z d z w b V m l / 3 5 0 D r J C h I D

1 < 2 <
6 e w B O / r F 7 s k o Y 4 = < / l a t e x i t > f a < l a t e x i t s h a 1 _ b a s e 6 4 = " O X X R E O K l f e l D g R k 9 C 3 K S V c m N / x s = " > A A A C E H i c b V D L S g M x F L 2 p r 1 p f V Z d u g k V w V W Z E 0 G X B j c u K 9 g H t U D J p p g 1 N M k O S E c r Q T 3 D j Q n / F n b j 1 D / w T l 2 b a W d j q g c D h n H u 5 J y d M B D f W 8 7 5 Q a W 1 9 Y 3 O r v F 3 Z 2 d 3 b P 6 g e H r V N n G r K W j Q W s e 6 G x D D B F W t Z b g X r J p o R G Q r W C S c 3 u d 9 5 Z N r w W D 3 Y a c I C S U a K R 5 w S 6 6 T 7 a E A G 1 Z p X 9 + b A f 4 l f k B o U a A 6 q 3 / 1 h T F P J l K W C G N P z v c Q G G d G W U 8 F m l X 5 q W E L o h I x Y z 1 F F J D N B N o 8 6 w 2 d O G e I o 1 u 4 p i + f q 7 4 2 M S G O m M n S T k t i x W f V y 8 V 8 v V 7 S J z N L 9 L J Q r e W x 0 H W R c J a l l i i 7 i R K n A N s Z 5 O 3 j I N a N W T B 0 h V H P 3 I 0 z H R B N q X Y c V V 5 W / W s x f 0 r 6 o + 4 7 f X d Y a 3 a K 0 M p z A K Z y D D 1 f Q g F t o Q g s o j O A J X u A V P a M 3 9 I 4 + F q M l V O w c w x L Q 5 w / d a Z 3 d < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 64 = " O X X R E O K l f e l D g R k 9 C 3 K S V c m N / x s = " > A A A C E H i c b V D L S g M x F L 2 p r 1 p f V Z d u g k V w V W Z E 0 G X B j c u K 9 g H t U D J p p g 1 N M k O S E c r Q T 3 D j Q n / F n b j 1 D / w T l 2 b a W d j q g c D h n H u 5 J y d M B D f W 8 7 5 Q a W 1 9 Y 3 O r v F 3 Z 2 d 3 b P 6 g e H r V N n G r K W j Q W s e 6 G x D D B F W t Z b g X r J p o R G Q r W C S c 3 u d 9 5 Z N r w W D 3 Y a c I C S U a K R 5 w S 6 6 T 7 a E A G 1 Z p X 9 + b A f 4 l f k B o U a A 6 q 3 / 1 h T F P J l K W C G N P z v c Q G G d G W U 8 F m l X 5 q W E L o h I x Y z 1 F F J D N B N o 8 6 w 2 d O G e I o 1 u 4 p i + f q 7 4 2 M S G O m M n S T k t i x W f V y 8 V 8 v V 7 S J z N L 9 L J Q r e W x 0 H W R c J a l l i i 7 i R K n A N s Z 5 O 3 j I N a N W T B 0 h V H P 3 I 0 z H R B N q X Y c V V 5 W / W s x f 0 r 6 o + 4 7 f X d Y a 3 a K 0 M p z A K Z y D D 1 f Q g F t o Q g s o j O A J X u A V P a M 3 9 I 4 + F q M l V O w c w x L Q 5 w / d a Z 3 d < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O X X R E O K l f e l D g R k 9 C 3 K S V c m N / x s = " > A A A C E H i c b V D L S g M x F L 2 p r 1 p f V Z d u g k V w V W Z E 0 G X B j c u K 9 g H t U D J p p g 1 N M k O S E c r Q T 3 D j Q n / F n b j 1 D / w T l 2 b a W d j q g c D h n H u 5 J y d M B D f W 8 7 5 Q a W 1 9 Y 3 O r v F 3 Z 2 d 3 b P 6 g e H r V N n G r K W j Q W s e 6 G x D D B F W t Z b g X r J p o R G Q r W C S c 3 u d 9 5 Z N r w W D 3 Y a c I C S U a K R 5 w S 6 6 T 7 a E A G 1 Z p X 9 + b A f 4 l f k B o U a A 6 q 3 / 1 h T F P J l K W C G N P z v c Q G G d G W U 8 F m l X 5 q W E L o h I x Y z 1 F F J D N B N o 8 6 w 2 d O G e I o 1 u 4 p i + f q 7 4 2 M S G O m M n S T k t i x W f V y 8 V 8 v V 7 S J z N L 9 L J Q r e W x 0 H W R c J a l l i i 7 i R K n A N s Z 5 O 3 j I N a N W T B 0 h V H P 3 I 0 z H R B N q X Y c V V 5 W / W s x f 0 r 6 o + 4 7 f X d Y a 3 a K 0 M p z A K Z y D D 1 f Q g F t o Q g s o j O A J X u A V P a M 3 9 I 4 + F q M l V O w c w x L Q 5 w / d a Z 3 d < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O X X R E O K l f e l D g R k 9 C 3 K S V c m N / x s = " > A A A C E H i c b V D L S g M x F L 2 p r 1 p f V Z d u g k V w V W Z E 0 G X B j c u K 9 g H t U D J p p g 1 N M k O S E c r Q T 3 D j Q n / F n b j 1 D / w T l 2 b a W d j q g c D h n H u 5 J y d M B D f W 8 7 5 Q a W 1 9 Y 3 O r v F 3 Z 2 d 3 b P 6 g e H r V N n G r K W j Q W s e 6 G x D D B F W t Z b g X r J p o R G Q r W C S c 3 u d 9 5 Z N r w W D 3 Y a c I C S U a K R 5 w S 6 6 T 7 a E A G 1 Z p X 9 + b A f 4 l f k B o U a A 6 q 3 / 1 h T F P J l K W C G N P z v c Q G G d G W U 8 F m l X 5 q W E L o h I x Y z 1 F F J D N B N o 8 6 w 2 d O G e I o 1 u 4 p i + f q 7 4 2 M S G O m M n S T k t i x W f V y 8 V 8 v V 7 S J z N L 9 L J Q r e W x 0 H W R c J a l l i i 7 i R K n A N s Z 5 O 3 j I N a N W T B 0 h V H P 3 I 0 z H R B N q X Y c V V 5 W / W s x f 0 r 6 o + 4 7 f X d Y a 3 a K 0 M p z A K Z y D D 1 f Q g F t o Q g s o j O A J X u A V P a M 3 9 I 4 + F q M l V O w c w x L Q 5 w / d a Z 3 d < / l a t e x i t > f a < l a t e x i t s h a 1 _ b a s e 6 4 = " O X X R E O K l f e l D g R k 9 C 3 K S V c m N / x s = " > A A A C E H i c b V D L S g M x F L 2 p r 1 p f V Z d u g k V w V W Z E 0 G X B j c u K 9 g H t U D J p p g 1 N M k O S E c r Q T 3 D j Q n / F n b j 1 D / w T l 2 b a W d j q g c D h n H u 5 J y d M B D f W 8 7 5 Q a W 1 9 Y 3 O r v F 3 Z 2 d 3 b P 6 g e H r V N n G r K W j Q W s e 6 G x D D B F W t Z b g X r J p o R G Q r W C S c 3 u d 9 5 Z N r w W D 3 Y a c I C S U a K R 5 w S 6 6 T 7 a E A G 1 Z p X 9 + b A f 4 l f k B o U a A 6 q 3 / 1 h T F P J l K W C G N P z v c Q G G d G W U 8 F m l X 5 q W E L o h I x Y z 1 F F J D N B N o 8 6 w 2 d O G e I o 1 u 4 p i + f q 7 4 2 M S G O m M n S T k t i x W f V y 8 V 8 v V 7 S J z N L 9 L J Q r e W x 0 H W R c J a l l i i 7 i R K n A N s Z 5 O 3 j I N a N W T B 0 h V H P 3 I 0 z H R B N q X Y c V V 5 W / W s x f 0 r 6 o + 4 7 f X d Y a 3 a K 0 M p z A K Z y D D 1 f Q g F t o Q g s o j O A J X u A V P a M 3 9 I 4 + F q M l V O w c w x L Q 5 w / d a Z 3 d < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O X X R E O K l f e l D g R k 9 C 3 K S V c m N / x s = " > A A A C E H i c b V D L S g M x F L 2 p r 1 p f V Z d u g k V w V W Z E 0 G X B j c u K 9 g H t U D J p p g 1 N M k O S E c r Q T 3 D j Q n / F n b j 1 D / w T l 2 b a W d j q g c D h n H u 5 J y d M B D f W 8 7 5 Q a W 1 9 Y 3 O r v F 3 Z 2 d 3 b P 6 g e H r V N n G r K W j Q W s e 6 G x D D B F W t Z b g X r J p o R G Q r W C S c 3 u d 9 5 Z N r w W D 3 Y a c I C S U a K R 5 w S 6 6 T 7 a E A G 1 Z p X 9 + b A f 4 l f k B o U a A 6 q 3 / 1 h T F P J l K W C G N P z v c Q G G d G W U 8 F m l X 5 q W E L o h I x Y z 1 F F J D N B N o 8 6 w 2 d O G e I o 1 u 4 p i + f q 7 4 2 M S G O m M n S T k t i x W f V y 8 V 8 v V 7 S J z N L 9 L J Q r e W x 0 H W R c J a l l i i 7 i R K n A N s Z 5 O 3 j I N a N W T B 0 h V H P 3 I 0 z H R B N q X Y c V V 5 W / W s x f 0 r 6 o + 4 7 f X d Y a 3 a K 0 M p z A K Z y D D 1 f Q g F t o Q g s o j O A J X u A V P a M 3 9 I 4 + F q M l V O w c w x L Q 5 w / d a Z 3 d < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O X X R E O K l f e l D g R k 9 C 3 K S V c m N / x s = " > A A A C E H i c b V D L S g M x F L 2 p r 1 p f V Z d u g k V w V W Z E 0 G X B j c u K 9 g H t U D J p p g 1 N M k O S E c r Q T 3 D j Q n / F n b j 1 D / w T l 2 b a W d j q g c D h n H u 5 J y d M B D f W 8 7 5 Q a W 1 9 Y 3 O r v F 3 Z 2 d 3 b P 6 g e H r V N n G r K W j Q W s e 6 G x D D B F W t Z b g X r J p o R G Q r W C S c 3 u d 9 5 Z N r w W D 3 Y a c I C S U a K R 5 w S 6 6 T 7 a E A G 1 Z p X 9 + b A f 4 l f k B o U a A 6 q 3 / 1 h T F P J l K W C G N P z v c Q G G d G W U 8 F m l X 5 q W E L o h I x Y z 1 F F J D N B N o 8 6 w 2 d O G e I o 1 u 4 p i + f q 7 4 2 M S G O m M n S T k t i x W f V y 8 V 8 v V 7 S J z N L 9 L J Q r e W x 0 H W R c J a l l i i 7 i R K n A N s Z 5 O 3 j I N a N W T B 0 h V H P 3 I 0 z H R B N q X Y c V V 5 W / W s x f 0 r 6 o + 4 7 f X d Y a 3 a K 0 M p z A K Z y D D 1 f Q g F t o Q g s o j O A J X u A V P a M 3 9 I 4 + F q M l V O w c w x L Q 5 w / d a Z 3 d < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " O X X R E O K l f e l D g R k 9 C 3 K S V c m N / x s = " > A A A C E H i c b V D L S g M x F L 2 p r 1 p f V Z d u g k V w V W Z E 0 G X B j c u K 9 g H t U D J p p g 1 N M k O S E c r Q T 3 D j Q n / F n b j 1 D / w T l 2 b a W d j q g c D h n H u 5 J y d M B D f W 8 7 5 Q a W 1 9 Y 3 O r v F 3 Z 2 d 3 b P 6 g e H r V N n G r K W j Q W s e 6 G x D D B F W t Z b g X r J p o R G Q r W C S c 3 u d 9 5 Z N r w W D 3 Y a c I C S U a K R 5 w S 6 6 T 7 a E A G 1 Z p X 9 + b A f 4 l f k B o U a A 6 q 3 / 1 h T F P J l K W C G N P z v c Q G G d G W U 8 F m l X 5 q W E L o h I x Y z 1 F F J D N B N o 8 6 w 2 d O G e I o 1 u 4 p i + f q 7 4 2 M S G O m M n S T k t i x W f V y 8 V 8 v V 7 S J z N L 9 L J Q r e W x 0 H W R c J a l l i i 7 i R K n A N s Z 5 O 3 j I N a N W T B 0 h V H P 3 I 0 z H R B N q X Y c V V 5 W / W s x f 0 r 6 o + 4 7 f X d Y a 3 a K 0 M p z A K Z y D D 1 f Q g F t o Q g s o j O A J X u A V P a M 3 9 I 4 + F q M l V O w c w x L Q 5 w / d a Z 3 d < / l a t e x i t > a l a t e x i t s h a 1 _ b a s e 6 4 = " L B J 0 b 2 d O m F t U X U E l 1 B q u A W l K e / I = " > A A A C E H i c b Z D L S g M x G I X / q b d a b 1 W X b o J F c F V m R N B l w Y 3 L i v Y C 7 V A y a a Y N z W V I M k I Z + g h u X O i r u F O 3 v o F v 4 t J M 2 4 V t / S H w c c 4 f c n K i h D N j f f / b K 6 y t b 2 x u F b d L O 7 t 7 + w f l w 6 O m U a k m t E E U V 7 o d Y U M 5 k 7 R h m e W 0 n W i K R c R p K x r d 5 H 7 r k W r D l H y w 4 4 S G A g 8 k i x n B 1 k n 3 u B f 0 y h W / 6 k 8 H r U I w h w r M p 9 4 r / 3 T 7 i q S C S k s 4 N q Y T + I k N M 6 w t I 5 x O S t 3 U 0 A S T E R 7 Q j k O J B T V h N o 0 6 Q W d O 6 a N Y a X e k R V P 1 7 4 0 M C 2 P G I n K b A t u h W f Z y 8 V 8 v V 7 S J z c L 7 W S S W 8 t j 4 O s y Y T F J L J Z n F i V O O r E J 5 O 6 j P N C W W j x 1 g o p n 7 E S J D r D G x r s O S q y p Y L m Y V m h f V w P H d Z a X W + p i V V o Q T O I V z C O A K a n A L d W g A g Q E 8 w Q u 8e s / e m / f u f c 5 W C 9 6 8 6 G N Y G O / r F 7 3 g n m c = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " L B J 0 b 2d O m F t U X U E l 1 B q u A W l K e / I = " > A A A C E H i c b Z D L S g M x G I X / q b d a b 1 W X b o J F c F V m R N B l w Y 3 L i v Y C 7 V A y a a Y N z W V I M k I Z + g h u X O i r u F O 3 v o F v 4 t J M 2 4 V t / SH w c c 4 f c n K i h D N j f f / b K 6 y t b 2 x u F b d L O 7 t 7 + w f l w 6 O m U a k m t E E U V 7 o d Y U M 5 k 7 R h m e W 0 n W i K R c R p K x r d 5 H 7 r k W r D l H y w 4 4 S G A g 8 k i x n B 1 k n 3 u B f 0 y h W / 6 k 8 H r U I w h w r M p 9 4 r / 3 T 7 i q S C S k s 4 N q Y T + I k N M 6 w t I 5 x O S t 3 U 0 A S T E R 7 Q j k O J B T V h N o 0 6 Q W d O 6 a N Y a X e k R V P 1 7 4 0 M C 2 P G I n K b A t u h W f Z y 8 V 8 v V 7 S J z c L 7 W S S W 8 t j 4 O s y Y T F J L J Z n F i V O O r E J 5 O 6 j P N C W W j x 1 g o p n 7 E S J D r D G x r s O S q y p Y L m Y V m h f V w P H d Z a X W + p i V V o Q T O I V z C O A K a n A L d W g A g Q E 8 w Q u 8 e s / e m / f u f c 5 W C 9 6 8 6 G N Y G O / r F 7 3 g n m c = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " L B J 0 b 2 d O m F t U X U E l 1 B q u A W l K e / I = " > A A A C E H i c b Z D L S g M x G I X / q b d a b 1 W X b o J F c F V m R N B l w Y 3 L i v Y C 7 V A y a a Y N z W V I M k I Z + g h u X O i r u F O 3 v o F v 4 t J M 2 4 V t / S H w c c 4 f c n K i h D N j f f / b K 6 y t b 2 x u F b d L O 7 t 7 + w f l w 6 O m U a k m t E E U V 7 o d Y U M 5 k 7 R h m e W 0 n W i K R c R p K x r d 5 H 7 r k W r D l H y w 4 4 S G A g 8 k i x n B 1 k n 3 u B f 0 y h W / 6 k 8 H r U I w h w r M p 9 4 r / 3 T 7 i q S C S k s 4 N q Y T + I k N M 6 w t I 5 x O S t 3 U 0 A S T E R 7 Q j k O J B T V h N o 0 6 Q W d O 6 a N Y a X e k R V P 1 7 4 0 M C 2 P G I n K b A t u h W f Z y 8 V 8 v V 7 S J z c L 7 W S S W 8 t j 4 O s y Y T F J L J Z n F i V O O r E J 5 O 6 j P N C W W j x 1 g o p n 7 E S J D r D G x r s O S q y p Y L m Y V m h f V w P H d Z a X W + p i V V o Q T O I V z C O A K a n A L d W g A g Q E 8 w Q u 8 e s / e m / f u f c 5 W C 9 6 8 6 G N Y G O / r F 7 3 g n m c = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " L B J 0 b 2 d O m F t U X U E l 1 B q u A W l K e / I = " > A A A C E H i c b Z D L S g M x G I X / q b d a b 1 W X b o J F c F V m R N B l w Y 3 L i v Y C 7 V A y a a Y N z W V I M k I Z + g h u X O i r u F O 3 v o F v 4 t J M 2 4 V t / S H w c c 4 f c n K i h D N j f f / b K 6 y t b 2 x u F b d L O 7 t 7 + w f l w 6 O m U a k m t E E U V 7 o d Y U M 5 k 7 R h m e W 0 n W i K R c R p K x r d 5 H 7 r k W r D l H y w 4 4 S G A g 8 k i x n B 1 k n 3 u B f 0 y h W / 6 k 8 H r U I w h w r M p 9 4 r / 3 T 7 i q S C S k s 4 N q Y T + I k N M 6 w t I 5 x O S t 3 U 0 A S T E R 7 Q j k O J B T V h N o 0 6 Q W d O 6 a N Y a X e k R V P 1 7 4 0 M C 2 P G I n K b A t u h W f Z y 8 V 8 v V 7 S J z c L 7 W S S W 8 t j 4 O s y Y T F J L J Z n F i V O O r E J 5 O 6 j P N C W W j x 1 g o p n 7 E S J D r D G x r s O S q y p Y L m Y V m h f V w P H d Z a X W + p i V V o Q T O I V z C O A K a n A L d W g A g Q E 8 w Q u 8 e s / e m / f u f c 5 W C 9 6 8 6 G N Y G O / r F 7 3 g n m c = < / l a t e x i t > a l a t e x i t s h a 1 _ b a s e 6 4 = " z f C w h 4 e H 8 o O W Q S 9 E A 7 U f d Q X z k R A = " > A A A C E H i c b Z D L S g M x G I X / 8 V r r r e r S T b A I r s p M E X R Z c O O y o r 1 A O 5 R M m m l D c x m S j F C G P o I b F / o q 7 t S t b + C b u D T T d m F b f w h 8 n P O H n J w o 4 c x Y 3 / / 2 1 t Y 3 N r e 2 C z v F 3 b 3 9 g v z X v 3 P m e r a 9 6 8 6 B N Y G O / r F 7 + K n m g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " z f C w h 4 e H 8 o O W Q S 9 E A 7 U f d Q X z k R A = " > A A A C E H i c b Z D L S g M x G I X / 8 V r r r e r S T b A I r s p M E X R Z c O O y o r 1 A O 5 R M m m l D c x m S j F C G P o I b F / o q 7 t S t b + C b u D T T d m F b f w h 8 n P O H n J w o 4 c x Y 3 / / 2 1 t Y 3 N r e 2 C z v F 3 b 3 9 g v z X v 3 P m e r a 9 6 8 6 B N Y G O / r F 7 + K n m g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " z f C w h 4 e H 8 o O W Q S 9 E A 7 U f d Q X z k R A = " > A A A C E H i c b Z D L S g M x G I X / 8 V r r r e r S T b A I r s p M E X R Z c O O y o r 1 A O 5 R M m m l D c x m S j F C G P o I b F / o q 7 t S t b + C b u D T T d m F b f w h 8 n P O H n J w o 4 c x Y 3 / / 2 1 t Y 3 N r e 2 C z v F 3 b 3 9 g v z X v 3 P m e r a 9 6 8 6 B N Y G O / r F 7 + K n m g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " z f C w h 4 e H 8 o O W Q S 9 E A 7 U f d Q X z k R A = " > A A A C E H i c b Z D L S g M x G I X / 8 V r r r e r S T b A I r s p M E X R Z c O O y o r 1 A O 5 R M m m l D c x m S j F C G P o I b F / o q 7 t S t b + C b u D T T d m F b f w h 8 n P O H n J w o 4 c x Y 3 / / 2 1 t Y 3 N r e 2 C z v F 3 b 3 9 g

2 <
e s / e m / f u f c 5 W C 9 6 8 6 G N Y G O / r F 7 3 g n m c = < / l a t e x i t > a l a t e x i t s h a 1 _ b a s e 6 4 = " z f C w h 4 e H 8 o O W Q S 9 E A 7 U f d Q X z k R A = " > A A A C E H i c b Z D L S g M x G I X / 8 V r r r e r S T b A I r s p M E X R Z c O O y o r 1 A O 5 R M m m l D c x m S j F C G P o I b F / o q 7 t S t b + C b u D T T d m F b f w h 8 n P O H n J w o 4 c x Y 3 / / 2 1 t Y 3 N r e 2 C z v F 3 b 3 9 g v z X v 3 P m e r a 9 6 8 6 B N Y G O / r F 7 + K n m g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " z f C w h 4 e H 8 o O W Q S 9 E A 7 U f d Q X z k R A = " > A A A C E H i c b Z D L S g M x G I X / 8 V r r r e r S T b A I r s p M E X R Z c O O y o r 1 A O 5 R M m m l D c x m S j F C G P o I b F / o q 7 t S t b + C b u D T T d m F b f w h 8 n P O H n J w o 4 c x Y 3 / / 2 1 t Y 3 N r e 2 C z v F 3 b 3 9 g v z X v 3 P m e r a 9 6 8 6 B N Y G O / r F 7 + K n m g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " z f C w h 4 e H 8 o O W Q S 9 E A 7 U f d Q X z k R A = " > A A A C E H i c b Z D L S g M x G I X / 8 V r r r e r S T b A I r s p M E X R Z c O O y o r 1 A O 5 R M m m l D c x m S j F C G P o I b F / o q 7 t S t b + C b u D T T d m F b f w h 8 n P O H n J w o 4 c x Y 3 / / 2 1 t Y 3 N r e 2 C z v F 3 b 3 9 g v z X v 3 P m e r a 9 6 8 6 B N Y G O / r F 7 + K n m g = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " z f C w h 4 e H 8 o O W Q S 9 E A 7 U f d Q X z k R A = " > A A A C E H i c b Z D L S g M x G I X / 8 V r r r e r S T b A I r s p M E X R Z c O O y o r 1 A O 5 R M m m l D c x m S j F C G P o I b F / o q 7 t S t b + C b u D T T d m F b f w h 8 n P O H n J w o 4 c x Y 3 / / 2 1 t Y 3 N r e 2 C z v F 3 b 3 9 g v z X v 3 P m e r a 9 6 8 6 B N Y G O / r F 7 + K n m g = < / l a t e x i t > f v < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n z 7 n w 4 n 4 v R k l P s n M I S n K 9 f A G q d 8 g = = < / l a t e x i t > t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " > A A A C E H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u C G 5 c V n b b Q D i W T Z t r Q J D M k m U I Z + g l u X O i v u B O 3 / o F / 4 t J M O w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n 6 2 r u m f 5 w 3 W t 0 S l K K 8 M Z n M M l e H A D D b i H J v h A Y A j P 8 A p v z o v z 7 n w 4 n 4 v R k l P s n M I S n K 9 f A G q d 8 g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " > A A A C E H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u C G 5 c V n b b Q D i W T Z t r Q J D M k m U I Z + g l u X O i v u B O 3 / o F / 4 t J M O w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n 6 2 r u m f 5 w 3 W t 0 S l K K 8 M Z n M M l e H A D D b i H J v h A Y A j P 8 A p v z o v z 7 n w 4 n 4 v R k l P s n M I S n K 9 f A G q d 8 g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " > A A A C E H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u C G 5 c V n b b Q D i W T Z t r Q J D M k m U I Z + g l u X O i v u B O 3 / o F / 4 t J M O w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n 6 2 r u m f 5 w 3 W t 0 S l K K 8 M Z n M M l e H A D D b i H J v h A Y A j P 8 A p v z o v z 7 n w 4 n 4 v R k l P s n M I S n K 9 f A G q d 8 g = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " d 9 3 g H A 1 f W j e D S w t f J i F s p 7 a e 6 a 0 = " > A A A C E H i c b V D L S g M x F L 1 T X 7 W + q i 7 d B I v g q s y I o M u C G 5 c V n b b Q D i W T Z t r Q J D M k m U I Z + g l u X O i v u B O 3 / o F / 4 t J M O w v b e i B w O O d e 7 s k J E 8 6 0 c d 1 v p 7 S x u b W 9 U 9 6 t 7 O 0 f H B 5 V j 0 9 a O k 4 V o T 6 J e a w 6 I d a U M 0 l 9 w w y n n FIG. 1.A cubical element of the mass-spring model (MSM).(a) The mass of the cube m is distributed equally among the 8 cubes that meet at one corner.The edge and diagonal springs have a stiffness of k/4 and k/2, respectively.(b) A Kelvin-Voigt type damping system is used where the edge and diagonal dampers have damping constant of c/4 and c/2, respectively.(c) To add volume conservation to the model, a penalty pressure is applied to the faces of the cube.This pressure results in forces at each corner.(d) The compressive active contraction is applied along the fiber direction in each cube and based on the distances a1 and a2 with bilinear interpolation, it has been redistributed on the masses.

A
< l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8= " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g Fe 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / 8 w y u 8 O Q / O i / P u f M x b C 8 4 i w m P 4 A + f z B 7 0 T j V M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8 = " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / 8 w y u 8 O Q / O i / P u f M x b C 8 4 i w m P 4 A + f z B 7 0 T j V M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8 = " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / 8 w y u 8 O Q / O i / P u f M x b C 8 4 i w m P 4 A + f z B 7 0 T j V M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8 = " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / 8 w y u 8 O Q / O i / P u f M x b C 8 4 i w m P 4 A + f z B 7 0 T j V M = < / l a t e x i t > tr(✏) < l a t e x i t s h a 1 _ b a s e 6 4 = " R M 7 U z A + M 4 S 2 4 K T M L Z H 6 q n 1 m Z N Q g = " > A A A B 8 3 i c b Z B L S w M x F I U z P m t 9 V V 2 6 C R a h b s q M C L o s u H F Z w T 6 g M 5 R M e t u G Z p I h u S O U o e C v c O N C E b f + G X f + G 9 P H Q l s P B D 7 O u S E 3 J 0 6 l s O j 7 3 9 7 a + s b m 1 n Z h p 7 i 7 t 3 9 w W D o 6 b l q d 5 d s C 4 E W 5 X y o f M M I 6 u p q I r I V j + 8 i o 0 L 6 u B X w 3 u r 8 q 1 + t O 8 j g I 5 J W e k Q g J y T W r k j t R J g 3 C S k m f y S t 6 8 z H v x 3 r 2 P + e i a t 6 j w h P y R 9 / k D 8 I O S D Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " R M 7 U z A + M 4 S 2 4 K T M L Z H 6 q n 1 m Z N Q g = " > A A A B 8 3 i c b Z B L S w M x F I U z P m t 9 V V 2 6 C R a h b s q M C L o s u H F Z w T 6 g M 5 R M e t u G Z p I h u S O U o e C v c O N C E b f + G X f + G 9 P H Q l s P B D 7 O u S E 3 J 0 6 l s O j 7 3 9 7 a + s b m 1 n Z h p 7 i 7 t 3 9 w W D o 6 b l q d 5 d s C 4 E W 5 X y o f M M I 6 u p q I r I V j + 8 i o 0 L 6 u B X w 3 u r 8 q 1 + t O 8 j g I 5 J W e k Q g J y T W r k j t R J g 3 C S k m f y S t 6 8 z H v x 3 r 2 P + e i a t 6 j w h P y R 9 / k D 8 I O S D Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " R M 7 U z A + M 4 S 2 4 K T M L Z H 6 q n 1 m Z N Q g = " > A A A B 8 3 i c b Z B L S w M x F I U z P m t 9 V V 2 6 C R a h b s q M C L o s u H F Z w T 6 g M 5 R M e t u G Z p I h u S O U o e C v c O N C E b f + G X f + G 9 P H Q l s P B D 7 O u S E 3 J 0 6 l s O j 7 3 9 7 a + s b m 1 n Z h p 7 i 7 t 3 9 w W D o 6 b l q d 5 d s C 4 E W 5 X y o f M M I 6 u p q I r I V j + 8 i o 0 L 6 u B X w 3 u r 8 q 1 + t O 8 j g I 5 J W e k Q g J y T W r k j t R J g 3 C S k m f y S t 6 8 z H v x 3 r 2 P + e i a t 6 j w h P y R 9 / k D 8 I O S D Q = = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " R M 7 U z A + M 4 S 2 4 K T M L Z H 6 q n 1 m Z N Q g = " > A A A B 8 3 i c b Z B L S w M x F I U z P m t 9 V V 2 6 C R a h b s q M C L o s u H F Z w T 6 g M 5 R M e t u G Z p I h u S O U o e C v c O N C E b f + G X f + G 9 P H Q l s P B D 7 O u S E 3 J 0 6 l s O j 7 3 9 7 a + s b m 1 n Z h p 7 i 7 t 3 9 w W D o 6 b l q d

FIG. 2 .
FIG. 2. A clockwise rotating spiral of active contraction and strain in a 2D isotropic elastic tissue.(a) Analytically derived active contraction field Ta(x, y).(b) The corresponding strain field Tr( )(x, y) is isotropic and closely follows the contraction field.(c) The mechanical phase map φ( )(x, y) reveals a mechanical phase singularity at the center of the rotor (red circle).(d) The complex amplitude A(x, y) vanishes at the position of the mechanical phase singularity.

A
< l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8= " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g Fe 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / / L D 5 2 S M 6 s M S B g r W 9 K Q X P 0 9 k d F I 6 0 k U 2 M 6 I m p F e 9 m b i f 1 4 3 N e G 1 n 3G Z p A Y l m y 8 K U 0 F M T G Z f k w F X y I y Y W E K Z 4 v Z W w k Z U U W Z s N i U b g r f 8 8 i p p X V Q 9 t + o 1 L i u 1 + t M 8 j i K c w C m c g w d X U I M 7 q E M T G C A8 w y u 8 O Q / O i / P u f M x b C 8 4 i w m P 4 A + f z B 7 0 T j V M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8 = " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / / L D 5 2 S M 6 s M S B g r W 9 K Q X P 0 9 k d F I 6 0 k U 2 M 6 I m p F e 9 m b i f 1 4 3 N e G 1 n 3G Z p A Y l m y 8 K U 0 F M T G Z f k w F X y I y Y W E K Z 4 v Z W w k Z U U W Z s N i U b g r f 8 8 i p p X V Q 9 t + o 1 L i u 1 + t M 8 j i K c w C m c g w d X U I M 7 q E M T G C A8 w y u 8 O Q / O i / P u f M x b C 8 4 i w m P 4 A + f z B 7 0 T j V M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8 = " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / / L D 5 2 S M 6 s M S B g r W 9 K Q X P 0 9 k d F I 6 0 k U 2 M 6 I m p F e 9 m b i f 1 4 3 N e G 1 n 3G Z p A Y l m y 8 K U 0 F M T G Z f k w F X y I y Y W E K Z 4 v Z W w k Z U U W Z s N i U b g r f 8 8 i p p X V Q 9 t + o 1 L i u 1 + t M 8 j i K c w C m c g w d X U I M 7 q E M T G C A8 w y u 8 O Q / O i / P u f M x b C 8 4 i w m P 4 A + f z B 7 0 T j V M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8 = " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b /

FIG. 3 .
FIG. 3. Focal active contraction and strain pattern in a 2D isotropic elastic tissue.(a) Circular ring-shaped focal pattern of active contraction Ta(x, y).(b) The corresponding strain field Tr( )(x, y) is isotropic and closely follows the circular active contraction field.(c) The mechanical phase map φ( )(x, y) exhibits a focal pattern and does not exhibit a phase singularity at any point.(d) The complex amplitude A(x, y) does not vanish at any point.

FIG. 4 .
FIG. 4. Strain spiral wave pattern in an anisotropic 2D medium with fibers aligned along the x-axis.(a) Clockwise rotating spiral wave pattern of active contraction Ta(x, y).(b) Strain along the fibers xx vanishes in the horizontally aligned arms.Red circles and black stars with white edge indicate the location of mechanical phase singularities calculated via the Hilbert transform and as the intersection of the xx = c contour (thin red dashed lines) and ∂t xx = 0 contours (thin blue lines), respectively, where c is the average of the minimum and maximum values of xx.(c) Phase map of the strain φ( xx(x, y)) shows one co-localized mechanical phase singularity at the center of the rotor pattern in (a) and two additional pairs in the top and bottom horizontally aligned spiral arms.(d) The amplitude A of the complex transformed signal vanishes in the vicinity to all mechanical singularities (black regions).(e) Time-series of xx for a points 1 and 2 in (b).(f) Phase diagram at points 1 and 2 from panel (e) derived via the Hilbert transform, where the horizontal axis is the original signal and the vertical axis is the signal shifted by π/2.

FIG. 5 .
FIG. 5. Focal wave in anisotropic medium with fibers aligned along x-axis.(a) Active contraction field Ta(x, y) of an elliptical ring-shaped focal pattern.(b) Strain xx is stronger in the vertical than the horizontal arms due to fiber anisotropy.Red circles and blue stars indicate the location of phase singularities calculated by the Hilbert and contour-intersection method as in Fig. 4, respectively.(c) The phase map shows mechanical singularities emerge in the absence of electrical singularities.(d) The signal amplitude A vanishes close to and around these singular points (black region).

FIG. 6 .
FIG.6.Mechanism of formation of unpaired mechanical singularities in anisotropic tissue with fibers aligned along the x-axis.Active tension (a) and strain xx (b) for a plane wave traveling in the +y-direction in a tissue with stress-free boundary conditions.Active tension (c) and strain xx (d) for a target wave expanding in a tissue with stress-free boundary conditions.Active tension (e) and strain xx (f) for a target wave expanding in a tissue with fixed boundary conditions (ux = uy = 0 on all boundaries).Red circles and blue stars indicate the location of phase singularities calculated by the Hilbert and contour-intersection method as in Figs.4 and 5, respectively.Panels (g) and (h) show the same strain patterns as (d) and (f), respectively, in deformed coordinates.

FIG. 7 .
FIG. 7. Number of mechanical phase singularities for a clockwise rotating spiral in a heart muscle with anisotropic fiber architecture (fibers are along x-axis) as a function of tissue size and tissue size to spiral wavelength ratio L/λs.Different vertical circles show the different number of phase singularities observed in one period.The intensity of the color in each circle exhibits the frequency of occurrence of that number of phase singularities.For instance, for (L/λs = 15) during one period 15, 17, and 19 phase singularities have been observed.The color density of the circles shows that most of the time there are 17 phase singularities in the domain.The red squares show the average number of singularities during one period.

FIG. 8 .
FIG. 8. Electromechanical phase singularities in the presence of a single scroll wave in a 3D tissue (1.2 × 1.2 × 1.2 cm).(a) All fibers are uniformly aligned along the x-axis.(b) Fibers orthotropically stacked and rotating by 90 • , parallel to xaxis at the bottom and parallel to y-axis at the top surface, respectively.Red lines represent lines of mechanical phase singularity or mechanical filaments and the black line at the center represents the electrical filament or scroll wave core.The iso-surface is presenting Ēf = Ēxx = 0 where Ēf is the normalized strain along the fiber that is used to calculate the phase.
t e x i t s h a 1 _ b a s e 6 4 = " B 9 w x I 5 O 7 / x p L a / t d L V 8 d 0 f k j / v g = " > A A A B 8 n i c b V D L S g N B E O y N r x h f U Y 9 e B o P g K e y K o B c h I t e x i t s h a 1 _ b a s e 6 4 = " D P N C 8 i v w 8 A y Y / 9 i n 6 f B n 4 0 n Y f 4 M = " > A A A B 6 n i c b V D L S g N B E O z x G e M r 6 t H L Y B A 8 h V 0 R 9 B j w 4 j F i X p A s o X c y m w y Z n V 1 m Z o W w B P w B L x 4 U 8 e o X e f N v n D w O m l j Q U F R 1 0 9 0 V p l I Y 6 3 n f Z G 1 9 Y 3 N r u 7 B T 3 N 3 b P z g s H R 0 3 T Z J p x h s s k Y l u h 2 i 4 F I o 3 r L C S t 1 P N M Q 4 l b 4 W j 2 6 n f e u T a i E T V 7 T j l Q Y w D J S L B 0 D r p o d 7 D X q n s V b w Z 6 C r x F 6 Q M C 9 R 6 p a 9 u P 2 F Z z J V l E o 3 p + F 5 q g x y 1 F U z y S b G b G Z 4 i G + G A d x x V G H M T 5 L N T J / T c K X 0 a J d q V s n S m / p 7 I M T Z m H I e u M 0 Y 7 N M v e V P z P 6 2 Q 2 u g l y o d L M c s X m i 6 J M U p v Q 6 d + 0 L z R n V o 4 d Q a a F u 5 W y I W p k 1 q V T d C H 4 y y + v k u Z l x f c q / v 1 V u V p 7 m s d R g F M 4 g w v w 4 R q q c A c 1 a A C D A T z D K 7 w R S V 7 I O / m Y t 6 6 R R Y Q n 8 A f k 8 w d H z o 4 6 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D P N C 8 i v w 8 A y Y / 9 i n 6 f B n 4 0 n Y f 4M = " > A A A B 6 n i c b V D L S g N B E O z x G e M r 6 t H L Y B A 8 h V 0 R 9 B j w 4 j F i X p A s o X c y m w y Z n V 1 m Z o W w B P w B L x 4 U 8 e o X e f N v n D w O m l j Q U F R 1 0 9 0 V p l I Y 6 3 n f Z G 1 9 Y 3 N r u 7 B T 3 N 3 b P z g s H R 0 3 T Z J p x h s s k Y l u h 2 i 4 F I o 3 r L C S t 1 P N M Q 4 l b 4 W j 2 6 n f e u T a i E T V 7 T j l Q Y w D J S L B 0 D r p o d 7 D X q n s V b w Z 6 C r x F 6 Q M C 9 R 6 p a 9 u P 2 F Z z J V l E o 3 p + F 5 q g x y 1 F U z y S b G b G Z 4 i G + G A d x x V G H M T 5 L N T J / T c K X 0 a J d q V s n S m / p 7 I M T Z m H I e u M 0 Y 7 N M v e V P z P 6 2 Q 2 u g l y o d L M c s X m i 6 J M U p v Q 6 d + 0 L z R n V o 4 d Q a a F u 5 W y I W p k 1 q V T d C H 4 y y + v k u Z l x f c q / v 1 V u V p 7 m s d R g F M 4 g w v w 4 R q q c A c1 a A C D A T z D K 7 w R S V 7 I O / m Y t 6 6 R R Y Q n 8 A f k 8 w d H z o 4 6 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D P N C 8 i v w 8 A y Y / 9 i n 6 f B n 4 0 n Y f 4 M = " > A A A B 6 n i c b V D L S g N B E O z x G e M r 6 t H L Y B A 8 h V 0 R 9 B j w 4 j F i X p A s o X c y m w y Z n V 1 m Z o W w B P w B L x 4 U 8 e o X e f N v n D w O m l j Q U F R 1 0 9 0 V p l I Y 6 3 n f Z G 1 9 Y 3 N r u 7 B T 3 N 3 b P z g s H R 0 3 T Z J p x h s s k Y l u h 2 i 4 F I o 3 r L C S t 1 P N M Q 4 l b 4 W j 2 6 n f e u T a i E T V 7 T j l Q Y w D J S L B 0 D r p o d 7 D X q n s V b w Z 6 C r x F 6 Q M C 9 R 6 p a 9 u P 2 F Z z J V l E o 3 p + F 5 q g x y 1 F U z y S b G b G Z 4 i G + G A d x x V G H M T 5 L N T J / T c K X 0 a J d q V s n S m / p 7 I M T Z m H I e u M 0 Y 7 N M v e V P z P 6 2 Q 2 u g l y o d L M c s X m i 6 J M U p v Q 6 d + 0 L z R n V o 4 d Q a a F u 5 W y I W p k 1 q V T d C H 4 y y + v k u Z l x f c q / v 1 V u V p 7 m s d R g F M 4 g w v w 4 R q q c A c 1 a A C D A T z D K 7 w R S V 7 I O / m Y t 6 6 R R Y Q n 8 A f k 8 w d H z o 4 6 < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " D P N C 8 i v w 8 A y Y / 9 i n 6 f B n 4 0 n Y f 4 M = " > A A A B 6 n i c b V D L S g N B E O z x G e M r 6 t H L Y B A 8 h V 0 R 9 B j w 4 j F i X p A s o X c y m w y Z n V 1 m Z o W w B P w B L x 4 U 8 e o X e f N v n D w O m l j Q U F R 1 0 9 0 V p l I Y 6 3 n f Z G 1 9 Y 3 N r u 7 B T 3 N 3 b P z g s H R 0 3 T Z J p x h s s k Y l u h 2 i 4 F I o 3 r L C S t 1 P N M Q 4 l b 4 W j 2 6 n f e u T a i E T V 7 T j l Q Y w D J S L B 0 D r p o d 7 D X q n s V b w Z 6 C r x F 6 Q M C 9 R 6 p a 9 u P 2 F Z z J V l E o 3 p + F 5 q g x y 1 F U z y S b G b G Z 4 i G + G A d x x V G H M T 5 L N T J / T c K X 0 a J d q V s n S m / p 7 I M T Z m H I e u M 0 Y 7 N M v e V P z P 6 2 Q 2 u g l y o d L M c s X m i 6 J M U p v Q 6 d + 0 L z R n V o 4 d Q a a F u 5 W y I W p k 1 q V T d C H 4 y y + v k u Z l x f c q / v 1 V u V p 7 m s d R g F M 4 g w v w 4 R q q c A c 1 a A C D A T z D K 7 w R S V 7 I O / m Y t 6 6 R R Y Q n 8 A f k 8 w d H z o 4 6 < / l a t e x i t > A < l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8 = " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / 8 w y u 8 O Q / O i / P u f M x b C 8 4 i w m P 4 A + f z B 7 0 T j V M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8 = " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / 8 w y u 8 O Q / O i / P u f M x b C 8 4 i w m P 4 A + f z B 7 0 T j V M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8 = " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b / 8 w y u 8 O Q / O i / P u f M x b C 8 4 i w m P 4 A + f z B 7 0 T j V M = < / l a t e x i t > < l a t e x i t s h a 1 _ b a s e 6 4 = " q M w F G f w t r n 0 1 y u z g 4 l B R Z a T 9 J u 8 = " > A A A B 6 H i c b V A 9 S w N B E J 2 L X z F + R S 1 t F o N g F e 5 E 0 D J i Y 5 m A + Y D k C H u b u W T N 3 t 6 x u y e E I 2 B v Y 6 G I r T / J z n / j 5 p J C E x 8 M P N 6 b Y W Z e k A i u j e t + O 4 W 1 9 Y 3 N r e J 2 a W d 3 b /

FIG. 10 .
FIG. 10.Electromechanical organization of filament-like phase singularities on the surface of a 3D slab with single scroll wave with small wavelength and fibers in x-direction.(a) Transmembrane voltage in dimensionless units.(b) Contraction field enslaved to voltage.(c) Strain depicted as largest compressive eigenstrain Exx.(d) Phase map φ(E) from the Hilbert transform of the normalized strain signal.(e) Amplitude map A(E) from Hilbert transform of the normalized strain signal.Amplitude vanishes in the neighbourhood of phase singular points.

FIG. 11 .
FIG. 11.Organization of electromechanical filaments in the presence of fiber rotation: (a-b) Single electrical scroll wave.(c-d)Composition of multiple scroll waves due to electrical wave break.(e-f) Focal wave pattern after application of a pulse at the center of the top surface.The isosurface is presenting Ēf = Ēxx = 0 where Ēf is the normalized strain along the fiber.Insets show co-localized electromechanical filaments.In the case that electrical waves remain stable one can see a clear co-localization between mechanical and electrical filaments.However, as scroll waves break, a more complex mechanism can be observed, but many of electrical filaments are still paired with a mechanical one.

FIG. 12 .
FIG. 12. Automated separation of paired (blue dots) from unpaired mechanical singularities (red crosses).Paired mechanical singularities co-localize with electrical singularities (black dots).The separation was performed by analyzing the temporal evolution of the strain surrounding a mechanical singularity.Phase singularities sampled over a short time interval (≈ 0.4s) on surface of 3D tissue during dynamics as shown in Figs.9-11: (a) Single scroll wave with fibers aligned along x-axis.(b) Single scroll wave with fibers rotated by 90 • (c) Scroll wave chaos with wave breakup.

FIG. 13 .
FIG. 13.Characteristic electrical impulse and contraction shapes exhibited by the electrophysiology model described in Sec.II B 1 and Appendix B. In this graph, a 1D representation of the equations is solved using the forward Euler method.In this simulation D = 1.1 cm 2 s −1 , Re = 0.8, dx = 3 × 10 −2 cm, and dt = 2.0 × 10 −4 s.Other parameters can be found in TableI.The left vertical axis corresponds to the variables U or v and the right vertical axis represents the dimensionless magnitude of the active contracting strain Ta.The horizontal axis is time in the units of milliseconds.

T = 1 s
→ L = 4.5 cm CT = 500 cm s × 1s = 500 cm T = 0.1 s → L = 4.5 cm CT = 500 cm s × 0.1s = 50 cm It is obvious that the assumption of L CT is valid for both cases.With ρ = 1000 kgm −3 the Young's modulus of the tissue is E = C 2 × ρ = 25 × 10 3 Nm −2 .Now, we need to confirm that with our assumptions, the size of the tissue is much less than the decay rate of the elastic waves.In our simulations the diffusivity of the tissue is D = 1.1 × 10 −4 m 2 s −1 .If we set D e in (33) to 1.1 × 10 −3 m 2 s −1 i.e., D e /D = 10 and further using the relations D e = η/ρ, δ = ηω/E, and ω = 2π/T , we obtain Eq. (47) 2.76 × 10 −4 has the unit of seconds.Substituting the excitation period over the range encompassing a normal heart rhythm and fibrillation, we obtain the estimates T = 1s → L = 4.5 cm CT πδ ≈ 2.89 × 10 6 cm T = 0.1s → L = 4.5 cm CT πδ ≈ 2.89 × 10 4 cm

τ el 2 F 6 D
substituting τ el = L 0 /C into (ρa 2 )/(Eτ 2 el ) will reduce Eq.(50) to:F = m∂ 2 t ū + c∂ t ū + kū (51)For the electrophysiology model, it is more convenient to adimensionalize time in Eq. (32) with τ = a 2 /D instead of τ el where D is the diffusion coefficient of membrane voltage in the electrophysiology model.In this case, the dimensionless elastic wave speed is τ /τ el and dimensionless damping constant is D e /D.The rescaled continuum elastodynamics equation becomes level of the MSM, Eq. (50) becomes τ normalizing time becomes important in setting the proper time step to carry out the numerical integration.In all the 3D simulations in this article we have set C = 5 ms −1 , dx = 3 × 10 −4 m, D = 1.1 × 10 −4 m 2 s −1 , and D e = 11 × 10 −4 m 2 s −1 .Those choices yield the time constants and adimensionalized elastic wave speed C τ = dx 2 D = (3 × 10 −4 ) 2 1.1 × 10 −4 = 8.18 × 10 −4 s .MECHANICAL PROPERTIES OF THE MSM

FIG. 14
FIG. 14. a long bar with aspect ratio 59/5 has been used to study the mechanical properties of the MSM.(a) Two equal and opposite displacements along the bar's long axis have been applied to the ends of the bar.This loading has been used to find the relation of the Young's Modulus and Poisson ratio with the penalty volumetric pressure and to find the stress-strain relation of the MSM.In this loading case, no boundary condition was introduced because at all times the bar is in force and moment equilibrium.(b) A periodic displacement has been applied to one end of the bar while the other end is kept fixed.This loading has been used to investigate the storage and loss modulus of the MSM.

FIG. 15 .
FIG. 15. (top)Young's modulus and (bottom) Poisson's ratio of the MSM at small strain = 0.1% and large strain = 15% as a function of dimensionless volumetric pressure p = p/E.Ēv is the dimensionless Young's Modulus which is equal to Ēv = Ev/E.Ev is from Eq. 11 and νv is from Eq. 12.

TABLE I .
Parameters for the electrophysiology and MSM models.n.u.stands for no unit.

TABLE II .
Parameters for MSM model used in Appendix D. n.u.stands for no unit.