Synchronization-based reconstruction of electromechanical wave dynamics in elastic excitable media

The beating of the heart is orchestrated by nonlinear waves of electrical excitation, which propagate rapidly through the heart muscle and trigger a release of intracellular calcium, which in turn fuels the contractions of cardiac muscle cells. During severe heart rhythm disorders such as atrial or ventricular fibrillation, the electrical excitation degenerates into complex, unorganized wave patterns, which cause highly irregular mechanical activations of the heart muscle. Within the ventricular wall, the onset of such arrhythmic electrical activity and progression into complex spatiotemporal electrical chaos is believed to be associated with three-dimensional electrical spiral or “scroll” vortex wave dynamics,^1–3 wave phenomena well known to be exhibited by excitable media.^4–6 However, to date, no imaging technique exists, which can provide adequate in-depth visualizations of the three-dimensional, intracardiac electrical wave phenomena, in particular, during cardiac fibrillation. As a result, the exact nature of the chaotic, fibrillatory electrical wave pattern remains poorly understood.

In the laboratory, the electrical waves can be imaged, directly and in great detail, on the heart’s surface using optical mapping.^2,7–10 Further, the epi- and endocardial electrical activity can be measured, at lower spatial resolutions, in human patients using electroanatomical catheter mapping¹¹ or contact electrode mapping.¹² The electrical activation can also be measured noninvasively, but rather indirectly, on the heart surface of patients using inverse electrocardiography.^13,14 However, inverse electrocardiographic techniques provide indirect measurements, the electrical activation patterns are calculated from a plurality of electrocardiograms measured on the body surface, and the technique’s accuracy and ability to resolve arrhythmias are under debate.^15–17 Transmural measurements of intracardiac electrical activity remain a scientific challenge, and it is debated whether the intramural activity can be inferred from observing the heart surface alone.^18,19 To date, the only way to measure the electrical activity directly within the heart muscle is to use plunge needle electrodes^20,21 or needle catheters.²² Panoramic visualizations of the intracardiac electrical waves can be obtained using transillumination imaging,^23,24 however, only as projections of the three-dimensional wave phenomena. In addition to direct measurement techniques, inverse computational approaches have been proposed to estimate the electrical activity within the cardiac muscle, from superficial or sparse electrical measurement data^25,26 as well as from mechanical deformation data.²⁷ 

In recent work, it was demonstrated that cardiac arrhythmias can be characterized by “mechanical” vortex wave phenomena, which can be imaged within the ventricular wall during tachycardia and fibrillation using high-speed 3D ultrasound.²⁸ Imaging the fibrillating, contracting ventricular wall simultaneously with electromechanical optical mapping^29,30 and high-speed 3D ultrasound, it was shown that action potential and calcium vortex waves, which appear on the contracting ventricular surface during fibrillation, lead to vortexlike deformations of the heart wall and, moreover, a filamentlike mechanical vortex structure, which appears to resemble the vortex filament structure of electrical scroll waves. The findings demonstrate that cardiac deformation mechanics can be highly correlated with the electrophysiological wave patterns even during arrhythmias and suggest that it may become possible to reconstruct intracardiac electrical activity from analyzing spatiotemporal mechanical wave patterns.

Here, we provide a numerical proof-of-principle that electrical excitation waves can be reconstructed throughout the volume of an elastic excitable medium by measuring its mechanical deformations. We employ a synchronization-based mechano-electrical data assimilation scheme, which integrates the mechanical measurement data into a numerical model of the measured elastic excitable medium. This electromechanical model then adapts to the mechanical input, such that its electrical part evolves as if it had initiated the mechanical contractions.

In previous work, it was shown that the dynamical activity evolving inside a spatially extended excitable medium can be measured and reconstructed in silico using a synchronization-based approach.³¹ Using a second excitable system that is a replication of the first excitable system, and establishing sensors, which couple the two systems by observing the dynamics in specific locations in the first system and inputting this information into the second system,³² it becomes possible to drive the second excitable system and allow it to evolve similarly to the first. In particular, the input via the controllers leads to electrical wave activity within the second system, which evolves in congruence with the electrical wave activity in the first, the two eventually becoming synchronized.

We adapted and modified this approach and developed a mechano-electrical data assimilation and reconstruction scheme by coupling two “elastic” excitable systems to each other; see Fig. 1. Instead of using the electrical part of the first system to drive the electrical part of the second system, we used the mechanical deformations exhibited by the elastic part of the first system to drive the electrical part of the second system, with the intent of achieving deformation in the second system that is similar to the deformation in the first. This was achieved by directly measuring the local elastic behavior in the first system and coupling the two systems, in the assumption that local contraction in the first system was an immediate response to electrical excitation. To confirm that both systems evolve similarly, we then compared their electrical excitation wave dynamics to each other. With this model-assisted approach, we show that it could be possible to drive and evolve a “virtual” numerical replication of the heart muscle tissue, within which electrical wave patterns, which cause elastic deformations of the tissue, can be reconstructed.

Excitable wave dynamics and elastic deformations of an anisotropic soft tissue were simulated in three-dimensional bulk-shaped simulation domains using a phenomenological electromechanical numerical model. To develop and test the reconstruction approach, two elastic excitable media, systems 1 and 2, were coupled to each other, where the first system created synthetic deformation data, which were then subsequently assimilated by the second system and used for the reconstruction. A schematic drawing of the reconstruction approach is shown in Fig. 1. The source code for the numerical model and reconstruction scheme is made available in a repository online (see the  Appendix).

Nonlinear waves of electrical excitation were simulated using the two-variable Aliev-Panfilov model³³ consisting of two coupled partial differential equations,

where $u$ and $v$ are normalized dimensionless excitation and recovery variables, respectively; $∇2u$ is a term describing the diffusive electrical coupling in between cells; $a$⁠, $b$⁠, and $k$ are model parameters; and

is a term modulating the restitution curve with the parameters $ϵ0$⁠, $μ1$⁠, and $μ2$⁠. The electrical model produces focal, planar, and reentrant spiral and scroll waves as well as spatiotemporal chaos. It was solved using forward Euler integration in a finite difference integration scheme. The parameter values used in this study for the electrical model and integration are given in Table I.

A coupling between electrics and mechanics was simulated using forward electromechanical coupling in the sense that excitation leads to mechanical contraction as previously described.³⁴ In each cell, excitation generates a temporary buildup of active stress $Ta$⁠. This stress generation was modeled using a third partial differential equation that was immediately coupled to the excitatory dynamic variable $u$ describing the evolution of the active stress,³⁴ 

Here, $kT$ is a parameter regulating the strength of the active stress buildup and $ϵ(u)$ models the twitchlike response,

Note that Eq. (5) is the corrected version⁴⁷ of the erratum in Eq. (23) of Nash and Panfilov.³⁴ As a result, nonlinear waves of electrical excitation are followed by waves of active contractile stress, where the active stress occurs along a fiber orientation (see Sec. II C). As the latter is a scalar-valued time-varying field $Ta(x,y,z,t)$⁠, equivalent to the excitatory and refractory dynamics in (1) and (2), the complete electrochemical simulation domain consists of three congruent simulation domains for excitatory, refractory, and active contractile stress dynamics.

Soft-tissue elasto-mechanics were simulated using a discrete particle system or mass-spring damper system, which allows the simulation of anisotropic contractions and deformations. The particles defining individual cells of the mechanical simulation domain are organized on a regular lattice forming hexahedral cells in a three-dimensional space. The barycenter $xb$ of each cell, which is defined as the average position of the cell’s 8 confining particles, is the origin of a set of 3 orthogonal spring pairs [see Fig. 2(b)]. To mimick cardiac muscle fiber anisotropy,³⁵ the spring orientations can be set individually to point into any direction $θ→(x,y,z)$ [see Fig. 2(b)]. Each spring intersects with two opposing faces of the hexahedral cell and is fixed to that position as the cell deforms. The intersection points $qj$ were computed in the undeformed configuration of the cell using a bilinear interpolation of the coordinates of the four vertices $x1$⁠, $x2$⁠, $x3$⁠, and $x4$ defining each face,

where $η$ and $ξ$ are bilinear coefficients. The intersection points were then used to compute the forces exerted by the springs onto the faces, and these forces were then redistributed via the four vertices to the neighboring cells and through the elastic system. The forces of each spring connecting the barycenter and an intersection point were computed as follows:

One of the three orthogonal spring pairs is set to be an active spring, which contracts due to a rise in the active stress variable $Ta$⁠,

Each cell of the mechanical model corresponds to one cell of the electrical model. Thus, active stress in one cell of the electrical model, see Eqs. (1)–(4), leads to a contraction of the active spring within one hexahedral cell of the elastic model. Additionally, structural springs were introduced along the edges between the vertices defining each cell [see also Fig. 2(b)],

In the above equations, $kj$⁠, $kf$⁠, and $kij$ are spring constants, $lj,0$ and $lf,0$ are rest lengths of the springs, and $cf$ is a scaling factor for the contraction strength (see Table I for parameter values).

Figure 2 shows a typical output of the electromechanical model. Electrical spiral wave activity (green), panel (a), induces a spiral-wave-like deformation of the elastic medium with areas of contraction (red) and dilation (blue) that align with the electrical wave fronts; see panel (c) and compare also with panel (d). Due to the presence of a muscle fiber architecture, see panel (b); contractions and deformations of the medium are stronger along the fiber orientation. For instance, in Figs. 2(a), 2(c)–2(e), and 3, the muscle fiber orientation is aligned uniformly, linearly transverse (see arrows) along the horizontal axis. In contrast, in Figs. 6 and 7, the medium retains a rotating fiber architecture. Due to the fiber anisotropy, the resulting deformation patterns are polarized, with maximal dilating and contracting strains occurring when waves travel along the fiber direction, while being diminished when waves travel in a perpendicular direction [see Figs. 2(c), 2(e), and 3].

To obtain elementary measures of deformation that could be analyzed and used within the reconstruction scheme, see Sec. II E, we performed a dimensionality reduction of tensorial kinematic observables of the mechanical configuration $χ1(t)$ of the elastic part of system 1, such that deformations of the elastic medium were expressed in terms of scalar-valued time-varying, three-dimensional deformation fields $E(x,t)$⁠. First, the right Cauchy-Green deformation tensor $C$ and the Green-Lagrangian deformation tensor $E$ were derived from the material prism given by the material coordinates defining each hexahedral cell of the elastic grid. Next, scalar-valued measures of deformation, such as volume or length changes, i.e., principal eigenvalues or other invariants of the deformation tensors, and their temporal derivatives were computed and used to describe deformation patterns throughout the medium. In this work, we used the volume change as a measure of deformation for the synchronization,

With this definition, negative values indicate local contraction and positive values indicate dilation of the hexahedral cells relative to their undeformed configuration. Throughout this study, contracted or contracting tissue is displayed in red, dilated or dilating tissue is displayed in blue, and undeformed tissue or tissue experiencing no deformation is displayed in white. As Fig. 2(c) shows, the deformation field can be displayed similarly in terms of volume change (left) or strain-rate (right). Qualitatively, both deformation patterns provide a similar picture of the deformation that was induced by a particular excitation pattern. Also, both patterns equally reflect the underlying fiber anisotropy.

To simulate experimental conditions, we added noise to the positions of the vertices $x~i=xi+ξ$ obtaining a noisy deformed mechanical configuration $χ~1(t)$ of system 1. Next, we computed the volume changes or strains and then smoothed the deformation fields using Gaussian smoothing [see Fig. 4(b)]. To counteract the polarization effects that occur due to the fiber architecture, see Figs. 2(c), 2(e), and 3, we further renormalized the resulting deformation time-series to the range $[−1,1]$ for each cell. This min-max normalization of a function $f(x,t)$ for a given point $x$ was accomplished through the scaling function,

where $mint(f(x))$ and $maxt(f(x))$ are the minimum and maximum of the time series $f(x,t)$⁠, respectively. As a result, the three-dimensional time-varying deformation field that was produced from the renormalized time-series contained a pronounced deformation field along all directions, including within the medium.

Two elastic excitable media, system 1 and system 2, were coupled using coupling elements consisting of sensors $Sijk$ in the elastic part of system 1 and controllers or drivers $Cijk$ in the electrical part of system 2 [see Fig. 1(b)]. The excitatory dynamics in systems 1 and 2 are described by

where $u1,2$ and $v1,2$ are the excitatory and refractory dynamic variables in the two systems, respectively; $F(⋅)$ denotes the right-hand side of the particular electrical model from Eq. (1); and $x~1$ refers to the obscured, noisy local kinematic quantities in system 1, such as positions and material lengths from which a particular measure of deformation was derived [see Eq. (10)]. The coupling signal $κ(⋅)$ enforces unidirectional coupling, as system 2 receives information from the sensors in system 1, but not vice versa. System 1 is the undisturbed dynamical system generating surrogate deformation data. By contrast, the evolution of system 2 is influenced by the coupling signal $κ(⋅)$⁠, which measures the locally averaged observed difference between the observed mechanical state in system 1 and the electrical state in system 2.

The spatial coupling scheme, illustrated in Fig. 1(b), places sensors $Sijk$ and controllers $Cijk$ sparsely throughout the medium of system 1 and system 2, respectively. Each coupling element consists of a pair of a block-shaped sensor $Sijk$ in system 1 and a block-shaped controller $Cijk$ with the same index $(i,j,k)$ at the same physical location in system 2 given in material coordinates. These sensors and controllers perform spatial averaging of the input (mechanical deformation of system 1) or control signal (excitation in system 2), respectively, on the voxels or simulation grid points within them. The sensors $Sijk$ analyze and process information about the elasto-mechanical state of system 1. For instance, the sensors could provide scalar-valued measures of deformation such as length or volume changes, i.e., principal eigenvalues or other invariants of the deformation tensor $E(x,y,z,t)$⁠, to the coupling element. The sensors and controllers are aligned on a regular lattice in material coordinates and are characterized by their size and the spacing in between them along each direction [see Fig. 1(b)]. The coupling signal is given by

where $κijk(t)$ is the coupling signal for all particles in a signal-controller pair and $u2$ is not directly manipulated in the space between controllers. This spatial coupling is an extension of a previously used coupling scheme for one-³² and two-dimensional³¹ spatially extended dynamical systems.

In addition to the sparse spatial alignment of coupling elements, the coupling is also modulated in time, illustrated in Fig. 1(b). The temporal coupling scheme aims to counteract (i) the delay between excitation and mechanical contraction and (ii) a sparse temporal sampling, which would be the case in real experimental imaging conditions. To correct the electromechanical delay, we employed the delay parameter $τ$ to shift the observations back in time. The coupling signal of a sensor-controller pair is based on the comparison of the averaged measure of the deformation within the sensor block in system 1 and the averaged excitation variable $u2$ within the controller,

Here, $s1(x,t)$ is an appropriate signal derived from a measure of the deformation of system 1, $N$ is the total number of cells or voxels contained in the sensor or controller, respectively, and $k±(t)$ is a time-dependent coupling strength. We make the assumption that the signal $s1(x,t)$ behaves in a similar way as the excitation variable $u1$⁠. This assumption holds when the excitation-contraction coupling mechanism is intact, electrical excitation is immediately converted into mechanical contraction, and a sequence of action potentials leads to a correlated and alternating sequence of dilation and contraction. The assumption is generally true for the normal heartbeat and is also motivated for ventricular tachycardia and fibrillation by recent experimental findings.²⁸ We chose $s1(x,t)=R(ΔV)(x,t)$⁠, the min-max normalized local volume change of each cell of the elastic medium [see Eqs. (10) and (11)].

The second part of the temporal coupling scheme is the choice of the coupling strength $k±(t)$ to compensate a sparse temporal resolution of the mechanical deformation [see Fig. 1(b)]. Anticipating the challenges of experimental limitations, the mechanical deformation of system 1 can only be observed by the synchronization every $n$th time step $tn,2n,3n,…$⁠. Figure 1(b) depicts the temporal coupling scheme and shows how the controllers use the input from system 1 at discrete and separate time points to adapt the electrical activity in system 2. In the reconstructions shown here, $n=20$ was chosen, and the coupling strength $k±(t)$ is modulated such that it decreases linearly with the temporal difference $Δt$ between the simulation time $t$ and the associated nearest observable data point $tobs$⁠,

Here, $k±$ are two coupling strength constants and $ks$ is a parameter to scale the linear decrease. The nearest observable data point $tobs$ is used for Eq. (15) to construct the continuous time series $s(x,t−τ)=s(x,tobs)$⁠; $tobs$ may also be in the future. For a given value $m∈N$⁠, the values of $Δt$ and $tobs$ between the two simulation time steps $tm⋅n+τ$ and $t(m+1)⋅n+τ$ are

Lastly, $k0±$ from Eq. (16) are two parameter constants; $k0+$ is used if $κijk(t)>0$ and $k0−$ if $κijk(t)<0$⁠, whereas $k±(Δt)$ is defined to be positive. This definition provides two coupling strength parameters, one if $u˙2$ is increased and the other one if $u˙2$ is decreased. This split may be beneficial the more the signal $s1(x,t)$ differs in behavior from $u1$⁠. For example, $k0−<k0+$ allows a more free-running synchronization, which counteracts regions where the signal $s1(x,t)$ is close to zero and $u1(x,t)$ is not [see Figs. 2(c), 2(e), and 3].

In summary, the reconstruction approach involves two coupled elastic excitable media: system 1 being the undisturbed system generating synthetic electrical and deformation data and system 2 coevolving with system 1 as it receives information about its mechanical state. By introducing sparse spatial and temporal observations, we intended to, firstly, mimick a real imaging situation, in which imaging speeds and spatial resolutions are constrained and, secondly, to take advantage of averaging effects and other intrinsic dynamical features of the model (see also Sec. IV). In terms of resolution, for instance, one rotation of a reentrant electromechanical wave could be imaged in real experimental imaging conditions with only 15–30 images with resolutions of only 20–100 pixels or voxels due to temporal or spatial resolution limits of the imaging device (see the  Appendix).

Elementary electrical excitation waves, such as planar, focal, and spiral waves, as well as three-dimensional electrical scroll wave dynamics and highly complex spatiotemporal electrical scroll wave chaos were reconstructed using the mechano-electric data assimilation approach [see Figs. 4–7]. In total, we performed more than 10 000 simulations and coupled elastic excitable system 2 to the dynamics that evolved in system 1 as described in Sec. II E. The simulations constitute systematic repetitions of four basic scenarios (thin spiral, thin focal wave, scroll wave with long wavelength, spatiotemporal chaos, Figs. 4–7) with different model parameters and varying synchronization parameters to optimize the reconstruction (see Fig. 10). The model and synchronization parameters used for the different scenarios are provided in Table I in the  Appendix.

Figure 4 illustrates how electrical spiral wave dynamics can be reconstructed from the analysis of mechanical spiral-wave-like deformations using the mechano-electric synchronization-based data assimilation approach. The upper row in panel (a) shows a single counterclockwise rotating electrical spiral wave in the electrical part of system 1. The medium is a thin three-dimensional bulk with two equally long sides (⁠$100×100×10$ cells), such that the dynamics are quasi-two-dimensional. The spiral’s wavelength is much smaller than the medium size, and it coils two times around its core region, which remains stationary at the medium center. The scenario is highly idealized and does not adequately describe wave dynamics during fibrillation in cardiac tissue. However, we use this example to illustrate the principle workings of the reconstruction technique. The lower row shows the response in the electrical part of system 2, as it assimilates the elasto-mechanics of system 1 using the sensors, controllers, and coupling as shown in Fig. 1(b). Within the reconstruction scheme, contractile deformation is perceived as the result of previous electrical excitation, and correspondingly, the controllers in system 2 initiate electrical excitation by pushing the excitatory dynamic variables in previous time steps above the excitation threshold. During the initial phase of the synchronization, the electrical activity corresponds to multiple wave sources that appear close to the mechanical wave front. Within a relatively short time, such wave sources merge and form electrical wave fronts, which become increasingly similar to the original wave pattern in system 1. Most importantly, due to the refractory dynamics of the excitable system, the waves start to propagate in the same direction as the original waves. As a result, contraction waves propagate through the medium and, in particular, the spatiotemporal evolution of the mechanical configuration of system 2 becomes increasingly similar to the one of system 1 [see panel (b)]. The times stated for each snapshot are normalized to one spiral rotation. In simulation time steps, the times correspond to $t=1190$⁠, $t=1688$⁠, $t=2450$⁠, and $t=3656$ from left to right. One spiral rotation corresponds to approximately $t=2950$ simulation time steps. The synchronization was enabled at $t=800$ and disabled at $t=10000$ simulation time steps. Within less than one and a half spiral rotations, the reconstructed electrical spiral wave in system 2 looks very similar to the original electrical spiral wave pattern in system 1. Also, the reconstructed elasto-mechanics that evolve in response to the electrical activation in the elastic part of system 2 are highly similar to the original elasto-mechanical pattern in system 1 [see panel (b)]. It is important to note that the reconstruction uses the obscured, noisy mechanical deformation pattern shown in panel (b) as input and was able to recover the original smooth deformation pattern. After the synchronization is turned off, the two spiral wave patterns look identical [see also Fig. 8(a)].

One critical observation, which is also further discussed in Secs. III E and III G, is that the synchronization procedure performs inefficiently when waves propagate along the perpendicular fiber direction. This inefficiency results from the muscle fiber orientation: the mechanical signal vanishes as the electrical waves propagate perpendicularly to the fiber orientation (see also Fig. 3). Subsequently, some parts of the spiral arms do not emerge at first in the electrics because the controllers do not inject sufficient excitation and the cells remain subthreshold. This behavior produces gaps in between sections of the electrical spiral arm [see $t=0.83$ in Fig. 4(a)]. Nevertheless, with time, the arms of the spiral wave grow toward each other and merge within less than two rotations because the dynamics in system 2 start to coevolve. The diffusivity of the excitable system takes over and establishes waves on its own, while the system receives simultaneously input from the controllers. The excitable system is, therefore, able to recover or correct missing or false information due to its intrinsic dynamics. In summary, it is possible to overcome complications that can be associated with anisotropy of the system, and one obtains two coevolving dynamics in the two elastic excitable systems. These two coevolving dynamics are surprisingly similar considering that the elasto-mechanics, which are driving the electrical dynamics, are substantially different from the original electrical wave dynamics. See also Sec. III G for a more detailed discussion of the residual mismatch between the two spirals and the reconstruction error.

The mechano-electric data assimilation approach can also recover focal wave patterns that originate from a point source (see Fig. 5). The electrical wave in system 1, see panel (a), originates from a single point at the center of a quasi two-dimensional flat bulk. The pulse is part of a sequence of pulses that are applied repetitively at the center of the medium. Each pulse is initiated after 5000 simulation time steps. The synchronization was started at $t=5000$ simulation time steps and continues throughout the sequence. As the tissue starts contracting at the center shortly after the application of the pulse, the sensors induce excitation in system 2 via the controllers [see panel (b)]. Consequently, the excitation develops quickly into a focus wave as well. This electrical focus wave receives further input from system 1 as it propagates outwards (⁠$t=0.32$⁠). Noise that leads to undesired electrical activity (⁠$t=0.78$⁠) at the medium boundary in system 2 is removed by subsequent pulses. System 2 solely reproduces the focal pattern from system 1 at times $t>1.0$⁠. The elastic bulk retains a linearly transverse fiber organization. Therefore, the mechanical contractile pattern is slightly anisotropic. As a result, the electrical focus wave in system 2 is slightly deformed in comparison with the one in system 1.

The mechano-electric reconstruction approach can furthermore also recover three-dimensional nonstationary scroll waves, which have a longer wavelength than the two previous examples. The electrical scroll wave shown in Fig. 6 retains a straight, slightly meandering-filament-like rotational core aligned along the $z$-axis located at the center of the bulk medium, see Fig. 6(c), and a wavelength that is in the same order as the medium size (⁠$λ≈d/3$⁠). The elastic bulk, which is thicker (⁠$Nz=25$ cells) than in the previous examples in Figs. 4 and 5, retains an orthotropic fiber organization with rotating sheets stacked along the $z$-axis. Due to the long wavelength, the deformations of the bulk are much larger and do not retain clear localized contraction wave patterns than in the examples in Figs. 4 and 5 [see Fig. 6(c)]. Nevertheless, even though there is no clear spiral or scroll-wave-like mechanical deformation pattern, the excitable system is able to recover the electrical scroll wave dynamics. It is also able to recover the electrical pattern in the presence of a complex fiber architecture. After 2–3 rotations, the two electrical scroll waves appear visually almost identical [see panel (a)].

Correspondingly, Fig. 6(d) shows the electrical vortex filaments of the original scroll wave in system 1 (gray) and the reconstructed scroll wave in system 2 (black). Both vortex filaments are colocalized beginning shortly after the start and throughout the synchronization procedure, and even after the synchronization was turned off, see also Fig. 6(e). Both vortex filaments are more or less straight, I-shaped, located at the center of the medium, and dissect the bulk in the $z$-direction. The original filament meanders slightly around a circular path (diameter $d≈15$ cells) and is perfectly straight, whereas the reconstructed filament meanders more along an elliptical and less regular path and bends and wraps at times around the original filament. The distance between the two filaments is small (⁠$4±2$ cells) and fluctuates only slightly [see Fig. 6(e)]. The reconstructed filament describes the core of the original scroll wave sufficiently accurate, given that the medium size is $100×100$ cells in the $x$- and $y$-direction and the scroll wave’s size corresponds approximately to the medium size. To characterize the mismatches between the two scroll waves more carefully, we computed their absolute difference as shown in Fig. 8(b); see also Sec. III G for a more detailed discussion of the reconstruction error. The times stated for the snapshots shown in Fig. 6 are normalized to one spiral rotation. In simulation time steps, the times in panel (b) correspond to $t=1320$⁠, $t=1610$⁠, $t=2100$⁠, and $t=4040$ from left to right and in panel (c) to $t=5505$⁠, $t=5925$⁠, and $t=7210$ from left to right. The synchronization was turned on at $t=1000$ and off at $t=10000$ simulation time steps. One rotation of the scroll wave corresponds to approximately $t=3100$ simulation time steps.

Even three-dimensional scroll wave chaos can be reconstructed sufficiently using the mechano-electric data assimilation approach (see Fig. 7). Panel (a) shows a comparison of the original electrical scroll wave dynamics that evolve in system 1 and the respective reconstructed electrical scroll wave dynamics that coevolve in system 2 as it receives input from system 1. Again, the three-dimensional bulk retains a rotating fiber architecture. Even though the original and reconstructed electrical wave patterns are not identical, they are qualitatively very similar; see also Movie S1 in the supplementary material. For instance, a double scroll wave or figure-of-eight pattern can be found in both systems at the top right center of the image. In system 2, the right arm of the reconstructed scroll wave has rotated slightly further than the original wave. However, aside from minor morphological differences, the two wave patterns have roughly the same shape, size, alignment, and topological charge. This is also demonstrated by the comparison of the original (gray) with the reconstructed (black) electrical vortex filaments in Figs. 7(c) and 7(d), which indicate the rotational core regions of the original and reconstructed scroll waves, respectively. The filaments appear in colocalized pairs (see pairs 1–3), which share the same vorticity and are in large parts coaligned. The distances between original and reconstructed filaments are small in comparison with the distances between the scroll waves, see the bottom right panel in (d), such that different scroll waves can clearly be distinguished from each other. However, despite the good agreement between the filaments, Fig. 8(c) shows that there are yet considerable mismatches between the two chaotic scroll wave dynamics when computing the absolute difference between the two spatial electrical patterns. Potential causes for these mismatches are discussed in Secs. III E and III G.

Overall, one can conclude that throughout the medium, the reconstruction is able to recover the essential parts of the scroll wave dynamics. In particular, the filaments demonstrate that the reconstruction recovers the three-dimensional topological organization of the electrical scroll wave dynamics. The reconstructed electrical wave dynamics are strikingly similar to the original electrical dynamics, given that they are produced from observations of obscured, noisy mechanical deformations [see Fig. 7(b)]. These are in principle fundamentally different from the electrical dynamics; see also Fig. 1(c). Due to the fiber anisotropy and nonlocality of the elasticity, the deformation pattern looks distorted or squashed when comparing it to the electrical wave pattern; see also Sec. III E. However, overall, the elastic deformations adumbrate the underlying electrical wave dynamics in that they exhibit similar length scales and the same direction of motion of the time-varying pattern, and this appears to be sufficient for the reconstruction.

The scroll wave chaos shown in Fig. 7 was preceded by a transient phase in which a single scroll wave emits focal and planar waves through the medium, which then develop into the chaotic dynamics through a cascade of wave breaks. In all cases, the reconstruction was able to recover the dynamics similarly as shown in Figs. 5 or 7; see also Movie S1 in the supplementary material.

Both fiber anisotropy and nonlocal elastic interactions produce significant alterations in the deformation patterns when comparing them with the electrical patterns. Time-series showing the elastic deformation of single tissue elements over time, see Figs. 1(c) and 3(b), do not exhibit the characteristic shapes of excitable waves with steep upstrokes, well-defined amplitudes, and diastolic intervals. The mechanical waves are correlated with the electrical wave sequences, see Fig. 1(c), but exhibit varying amplitudes and shapes, which depend on various factors, such as the spatial orientation of the waves with respect to the fiber orientation, their position within the medium, or contractile or tensile strains at a distance; see Figs. 1(c), 2(c) and 2(e), and 3(b). These phenomena pose a particular challenge to the reconstruction procedure. We introduced the renormalization, see Eq. (11) and Figs. 4(b) and 9, to compensate the amplitude variations of the mechanical waves. As discussed for the single spiral shown in Fig. 4, the renormalization helps to overcome problems associated with weakly deforming areas in the tissue that do not exhibit a notable sequence of dilation and contraction; see Figs. 2(e) and 3. By amplifying such faint signals, the topology of the original electrical wave pattern can be recovered: the diffusive electrical dynamics in system 2 close the blank areas and restore the wave fronts autonomously. As Fig. 7 demonstrates, this restoration mechanism works even with complicated and highly irregular scroll wave chaos.

To simulate experimental conditions, we added noise to the positions of the deformed mechanical lattice of system 1 [see Figs. 4(b), 6(c), 9, and 14]. In a real imaging situation, the noise could be evoked by noisy image data or the intrinsic noise of a tracking algorithm. As Figs. 4–7 demonstrate, the reconstruction technique performs robustly in the presence of noise and was able to recover elementary spiral and scroll waves as well as chaotic scroll wave dynamics. Even as the renormalization introduces additional noise and leads to substantial noise levels in areas with faint mechanical signals (the amplification of the faint signal also amplifies noise), the excitable part of system 2 is able to smooth out the noisy input from system 1 (see Fig. 9). However, Fig. 9 also shows that the renormalization of the faint signal may also cause a substantial phase shift between the original and reconstructed electrical signals [cf. Fig. 3(b)].

The performance of the reconstruction technique was measured and evaluated by comparing the electrical dynamics in two systems 1 and 2. To find optimal synchronization performance, we performed approximately 10 000 simulations and varied the different synchronization parameters, which would minimize the reconstruction error or mismatch of the dynamics between the two systems (see Table I in the  Appendix). In Figs. 6 and 7, we already demonstrated that electrical scroll vortex wave filaments can be reconstructed with sufficient accuracy by observing the mechanics. Here, we discuss the absolute mismatch between the original and reconstructed electrical patterns and the origins of systematic reconstruction errors.

Figure 8 shows the absolute difference (RMSE) between the original electrical dynamics in system 1 and the reconstructed electrical dynamics in system 2 during the different scenarios shown in Figs. 4–7. The error lies within the interval $[−1,1]$ and is displayed in violet and green, respectively (green displays remaining wave activity from system 1 and violet displays wave activity from system 2 when the two activities are subtracted from each other). Figure 8(a) shows the thin spiral in systems 1 and 2 after the synchronization was turned off (⁠$t=10000$ simulation time steps); see also Fig. 11. At this point, both spirals evolve independently from each other, and they look overall almost identical. They appear to rotate at the same angular speed and appear to be in phase. However, the two spirals are not perfectly congruent. The error map indicates slight offsets at the wave fronts and backs. It is important to point out that these mismatches particularly occur when the electrical waves propagate in a perpendicular direction to the fiber orientation. One can also see that the reconstructed wave slightly precedes the original wave (violet).

The residual mismatches in the electrics occur, even though the optimal synchronization parameters had already been chosen for the reconstructions shown in Figs. 4–7. For instance, we fine-tuned the synchronization parameter $τ$⁠, see Eq. (15), which controls the negative electromechanical delay. We introduced this parameter to compensate the electromechanical delay or delayed onset of contraction following the excitation in the model [see Fig. 1(b)]. Figures 10(a) and 10(b) show how the reconstruction error can be minimized for a particular negative electromechanical $τ$ for the thin spiral wave (⁠$τ=80$ simulation time steps or $0.03$ rotations) and thick scroll wave (⁠$τ=240$ simulation time steps or $0.08$ rotations), respectively. The error is stated as the mean root mean square error (mean RMSE). Selecting an optimal value for the negative electromechanical delay $τ$ during the mechano-electric reconstruction is essential as the reconstruction assumes that contractile deformation occurred due to previous excitation. Correspondingly, one has to introduce a coupling that modulates the electrical dynamics in system 2 in previous time steps to compensate the temporal delay that arises in both systems in between electrics and mechanics. If the reconstruction would not include a compensation of the electromechanical delay, the reconstructed dynamics would always slightly precede the original dynamics. Given the example in Fig. 4, the spiral in system 2 would have a slight uniform phase shift, and the error map would indicate a spiral wave pattern that corresponded to the mismatch between the two spirals. Therefore, by selecting an optimal negative electromechanical delay $τ$⁠, one can minimize a phase shift between systems 1 and 2 and reduce the global reconstruction error [see Figs. 10(a) and 10(b)].

However, even with an optimal negative electromechanical delay $τ$⁠, we could not fully compensate the reconstruction error at the wave fronts and backs shown in Fig. 8(a). The residual mismatches are an intrinsic feature of our reconstruction approach, and they are related to particular morphological features of the elastic wave pattern. More precisely, the weak mechanical signal that occurs in regions in which waves propagate in the perpendicular direction to the local fiber direction, see (4) in Fig. 3(b), is ambiguous and contains a phase shift with respect to the original electrical wave. This phase shift is consequently introduced in system 2 and causes the reconstruction errors at the wave fronts and backs. We conclude that such effects cannot be compensated by a single global parameter $τ$⁠.

To further optimize the synchronization, we varied the sensor spacing and sensor size, see Fig. 10(c), and found that for dynamics with smaller wavelengths, it is advantageous to use smaller sensors with smaller distances between the sensors, see Figs. 4 and 7, and correspondingly larger sensors and larger distances between the sensors for dynamics with larger wavelengths, see Fig. 6 and also Table I in the  Appendix. Figure 10(c) shows the parameter scan for the scroll wave with a large wavelength shown in Fig. 6. The minimal reconstruction error occurs along a linear basin of combinations of sensor sizes and sensor spacings with an absolute minimum at a sensor size of $3×3×2$ and a sensor spacing of $4×4×2$⁠. Within the linear basin, the minimal sensor spacing can scale slightly faster than the minimal sensor size. In comparison, the ideal sensor sizes and spacings were found to be $2×2×1$ and $1×1×1$⁠, respectively, for small wavelengths; see also Table I in the  Appendix.

Overall, the finite, nonvanishing reconstruction error shown in Fig. 8 highlights the discrepancies that remain between the original electrical and the reconstructed electrical dynamics despite the optimization of various synchronization parameters. Panels (b) and (c) show the reconstruction error for the single scroll wave shown in Fig. 6 and the scroll wave chaos shown in Fig. 7. The error does not vanish and is significant during the scroll wave chaos. With the single large scroll wave, the reconstruction error fluctuates heavily as the scroll wave rotates, see also Fig. 11(b), and it reaches levels in which there are marginal mismatches (bottom) and stronger mismatches (top), see Fig. 8(b), which can presumably be associated with boundary and finite size effects. During scroll wave chaos, see Fig. 8(c), the reconstruction error appears throughout the medium with positive and negative patches appearing close to wave fronts and backs. Such patches relate to the local deviations of the wavelets in the original and reconstructed system, respectively, as shown in Fig. 7(a).

In the previous examples, we performed a “twin experiment,” in which two systems 1 and 2 are identical with the same model equations and identical parameter sets. Nevertheless, to investigate whether our reconstruction technique can also be used in situations, in which the two systems are not completely identical, we coupled two systems with a different electrical parameter $a$⁠, which modulates the excitation threshold [see Eq. (1)]. Figure 12 demonstrates that the scroll wave shown in Fig. 6 can also be recovered even with a substantial parameter mismatch. We chose two parameter values, a higher (⁠$a=0.13$⁠, $Δa=0.03$⁠) and a slightly lower (⁠$a=0.0985$⁠, $Δa=−0.0015$⁠) parameter value $a2$ in system 2, while the original parameter $a1=0.1$ stayed the same in system 1. Figure 12(a) demonstrates that despite the large parameter mismatch (⁠$Δa=0.03$⁠), the original scroll wave in system 1 can be recovered in system 2 during continued coupling between systems 1 and 2 (⁠$t=2.83$⁠). Accordingly, in the upper panel of Fig. 12(b), it is shown that despite the large $Δa=0.03$⁠, both scroll waves are nearly congruent. The reconstructed scroll wave lags slightly behind the original scroll wave as indicated by the offsets at the wave’s fronts and backs (absolute difference or RMSE) [cf. Figs. 8 and 11(b)]. However, the lower panel in (b) shows that if the synchronization is turned off, the scroll wave in system 2 quickly develops (⁠$t=3.69$⁠) into a different scroll wave that is not following the dynamics in system 1 any longer. During the mechano-electric synchronization (⁠$t=2.83$⁠), the two electrical filaments are colocalized, cf. Fig. 6(d), whereas without synchronization (⁠$t=3.69$⁠), the two filaments no longer stay colocalized. The distance between the two filaments grows quickly over time, see Fig. 11(b) (dark-gray), as the reconstructed filament moves outwards toward the boundary of the bulk. After the coupling is turned off, the two scroll waves correspond after only one rotation (⁠$t=4.35$⁠) to two completely disparate scroll wave patterns. The quick dissociation between the two dynamics without coupling demonstrates that the parameter mismatch $Δa=a2−a1=0.03$ is substantial. In turn, panel (a) demonstrates that when the synchronization is on and system 2 receives continuous input from system 1, it can be kept coupled to system 1 even at substantial parameter mismatches. Figure 12(c) compares the congruency of the electrical scroll wave patterns after the synchronization was turned off for identical parameters $Δa=0$ (center) and the parameter mismatches $Δa=−0.0015$ (top) and $Δa=0.03$ (bottom). For $Δa=0$⁠, the absolute difference corresponds to the reconstruction error with identical systems 1 and 2. With small negative mismatch $Δa=−0.0015$⁠, the absolute difference is comparable as with $Δa=0$⁠, even though the reconstructed wave precedes the original wave slightly further than with matching parameters. With the larger positive mismatch, the reconstructed wave lags behind.

The data demonstrate that the mechano-electric synchronization can in principle recover the original scroll wave dynamics even with parameter mismatches between the two systems. This suggests that the reconstruction technique could potentially also work in a real imaging situation, in which the parameters of system 1 may be unknown and may have to be guessed in system 2.

In this study, we used a synchronization-based mechano-electric data assimilation approach to demonstrate that electrical excitation wave dynamics can be reconstructed within elastic excitable media by observing and analyzing the spatiotemporal mechanical deformation that the electrical excitation waves had caused. We provided a numerical proof-of-principle, which shows that it is possible to reconstruct various electrical excitation wave dynamics evolving in two- and three-dimensional anisotropic elastic excitable media, from elementary focal and planar waves to spiral and scroll waves, to even complicated spatiotemporal chaos, solely by observing soft-tissue mechanics. We used an approach that resembles the synchronization of two dynamical systems to each other,^31,32 where the dynamics in one system are measured by another system via a particular coupling between the two systems. The second system then adapts to the dynamics and begins to coevolve with the first system, to eventually create a replication of its full dynamical state. Our reconstruction technique can be compared to data assimilation approaches in weather prediction applications, in which observations of meteorological variables such as temperature, pressure, or wind speeds are used to initialize and evolve numerical weather forecast models.^36–39 In our adaptation of such data assimilation schemes, we coupled two electromechanical models, systems 1 and 2, and assumed that contractile deformation observed in the first system was caused in immediate response to previous electrical excitation. Correspondingly, excitation was initiated within the second system in locations, in which the first system contracted. As the excitation evolves in the second system, it initiates contractions and deformation, and soon the deformation of the second system begins to resemble the deformation of the first system. With continued coupling, the electrical and mechanical dynamics in the second system coevolve with the dynamics in the first system.

Given that our approach observes and assimilates solely spatiotemporal mechanical deformation patterns, it is remarkable that it is able to recover even complicated, highly irregular electrical scroll wave chaos (see Fig. 7). Comparing the two morphologies of the electrical and mechanical patterns, see Figs. 6 and 7, it is not always immediately evident, which electrical pattern underlies a particular mechanical pattern. Long-range elastic interactions, as well as finite size and boundary effects, create a complex relationship between electrics and mechanics. Even the single scroll wave in Fig. 6 does not simply induce a scroll-wave-like deformation pattern, but produces a superposition of local contractile strain and global deformation patterns. Figure 1(c) shows that the elastic waves exhibit, unlike the corresponding electrical waves, arbitrary shapes and amplitudes. Furthermore, they are not always strictly in phase with the electrical wave fronts, see also Fig. 9, and do not retain a refractory period. Lastly, the nonlocality of the elasticity and the fiber anisotropy produce mechanical patterns that appear to be highly distorted renditions of the electrical patterns. Nevertheless, the recovery of the electrical wave dynamics, even the complicated scroll wave chaos, is achieved as the second system intrinsically retains the same excitatory and refractory reaction-diffusion dynamics as the measured system. The mechanical observations “serve as a guiding ‘potential’ working with the dynamical rules”³⁹ of the electrical part of the second system to direct its overall electromechanical dynamics. The diffusive nature and local kinetics of the electrical dynamics introduce directionality and a refractory phase, within which further contractions are considered unphysiological and are prevented. As a result, it appears that it does not matter how exactly the excitation translates into a deformation pattern. As long as the characteristic length and time scales of the two wave dynamics are comparable, it is possible to recover an approximate reconstruction of the original electrical dynamics, in which the topological organization remains preserved.

Our reconstruction technique could be used to estimate and visualize intramural electrophysiological wave dynamics during the imaging of heart rhythm disorders, such as atrial fibrillation or ventricular tachycardia or other heart dysfunctions. The technique could be applied in combination with numerical motion tracking with high-speed 3D ultrasound, magnetic resonance imaging, or optical or other tomographic imaging techniques to measure the contractile motion and deformation of the heart at high speeds. The imaging and motion tracking data could then be assimilated in a numerical model, which dynamically adapts to the mechanical imaging data and simultaneously reconstructs the electrical activity that caused the contractile deformations. The numerical model would ideally correspond to an anatomically detailed virtual replication of the measured heart, be constructed dynamically in real-time from the images, and feature the heart’s fiber organization, electrophysiology, tissue mechanics as well as possible anatomical abnormalities or pathophysiologies, such as fibrotic or infarcted tissue regions.

Nowadays, with recent advancements in ultrasound technology, it is possible to image large parts of the ventricular walls in 3D at speeds (100–200 volumes per second) sufficient to resolve complex spatiotemporal mechano-electrical processes, such as three-dimensional mechanical vortex wave dynamics during arrhythmias.²⁸ At these imaging speeds, it is possible to obtain 15–30 volumetric images of one rapid arrhythmic contraction cycle (period of contraction 100–200 ms or dominant frequencies of 5–10 Hz during ventricular arrhythmias²⁸); see also the  Appendix. Also, modeling has advanced to a point such that large-scale cardiac simulations can now be performed in real-time on compact and moderately priced GPU-based computing architectures,⁴⁰ suggesting that it could become possible to combine sufficiently fast numerical models with our proposed imaging approach for real-time applications in the near future. Our data assimilation approach could also be combined with alternative, possibly complementary approaches for the estimation of excitation from cardiac muscle deformation²⁷ as well as with other purely electrical inverse computational techniques.^25,26

The central assumption in our study is that we expect that it will be feasible in future applications to employ a numerical model that “already” describes the dynamics of the observed system sufficiently accurate. Nevertheless, presently, it remains unclear as to whether computer models adequately describe the three-dimensional electrical wave phenomena during arrhythmias, as scroll wave dynamics have never been imaged in full inside the heart walls, and the models lack validation. Our method, however, could be employed to develop and validate computer models as the comparison of the model with the measurement data is part of the method’s principle workings. We believe that our approach will work independently from the choice of the particular numerical model and could potentially be performed with both phenomenological or with more realistic electrophysiological or mechanical models. In our study, we demonstrated that the mechano-electric reconstruction can be performed with two “detuned” systems with mismatching model parameters, see Sec. III H, suggesting that the technique is applicable in scenarios, in which the model parameters (system 2) are unknown and would have to be estimated and, as a result, could contain errors or differ from the true system parameters of the real observed system (heart).

A fundamental limitation of our reconstruction approach is associated with the biophysical principles of cardiac electromechanics: the reconstruction relies on the intact coupling of electrics and mechanics, and it would be compromised by a degeneration of the excitation-contraction coupling mechanism. To provide a proof-of-principle that mechanical deformation patterns can be used to retrieve electrical patterns, we employed elementary forward-coupling consisting of a direct conversion of electrical excitation into active contractile stress. This implementation is a crude simplification of cardiac electromechanics and not realistic. However, it served the purpose to provide a proof-of-principle and, more importantly, is also motivated by experimental evidence: it was shown that (i) intracellular calcium oscillations stay finite and correlated with transmembrane voltage oscillations during arrhythmias⁴¹ and (ii) that excitation-contraction coupling may remain largely intact during tachycardia and fibrillation²⁸ and (iii) produces electromechanical waves on the tissue scale,^28,42 which can in turn produce highly correlated spatiotemporal electromechanical dynamics even during ventricular fibrillation.²⁸ Our reconstruction technique may become compromised in the presence of electromechanical dissociation phenomena, i.e., a dissociation between voltage and calcium,^41,43–45 mechanical inhomogeneities, or other electromechanical dysfunctions,⁴⁶ and we aim to explore the effect of such phenomena onto the efficacy of our method in future research.

Further potential limitations of our study are associated with the specific implementation of the mechano-electric reconstruction in conjunction with our current model. We believe that, in principle, one could use many other mechanical observables or measures of deformation in the mechano-electric coupling scheme, i.e., principal eigenstrains, other invariants of the strain tensor, projections of length changes onto predefined orientations, angular changes, or also directly tensorial data itself, as long as one obtains a sequence of dilation and contraction that correlates with the sequence of electrical excitation, and as long as contraction can be uniquely identified, such that it triggers excitation at the correct times. Therefore, we anticipate that the reconstruction will work independently from the choice of a specific mechanical measure and will be successful as long as the overall deformation pattern that one constructs from the mechanical configuration is sufficiently similar to the electrical pattern. Figure 2 depicts that both the local volume change and strain-rate may similarly reflect the underlying electrical pattern. However, we used a mass-spring damper system in our study, which is not perfectly volume-preserving as it deforms. In compressible elastic media, the volume change is an appropriate measure to identify dilating and contracting regions. For instance, contraction and deformation waves in engineered cardiac tissue cultures may exhibit substantial density and volume changes. However, volume change is not suited for describing ventricular muscle tissue, which is considered to be incompressible. Nevertheless, one could define an appropriate mechanical signal that would adequately reflect the underlying excitation pattern. We aim to address these questions in a more systematic study in future research.

We provided a numerical proof-of-principle, which demonstrates that it is possible to reconstruct complex electrical wave dynamics, such as spiral and scroll waves in an elastic excitable cardiac tissue, solely by observing and analyzing spatiotemporal mechanical deformations. The reconstruction was enabled by integrating the data into a computer model of cardiac physiology with mechano-electric coupling. Our findings suggest that the approach could be used in combination with high-speed imaging, such as 3D ultrasound or magnetic resonance imaging, to noninvasively measure electrophysiological wave phenomena within the heart muscle, including the irregular activity during arrhythmias.

See the supplementary material for reconstruction of chaotic regime simulation shown in Figs. 1, 7, and 8 (Video S1); reconstruction of the thin spiral wave shown in Figs. 4 and 8 (Video S2); visualization of the simulation grid shown in Fig. 2(d) (Video S3); reconstruction of the thick spiral wave shown in Figs. 6 and 8 (Video S4); electrical filaments of the original and reconstruction system shown in Figs. 6(d) and 6(e) (Video S5); and electrical filaments of the original and reconstruction system shown in Figs. 7(c) and 7(d) (Video S6).

We would like to thank Ulrich Parlitz for fruitful discussions. The research was funded by the German Center for Cardiovascular Research (DZHK e.V.), partnersite Göttingen. We would like to thank the Max Planck Society for providing access to computing infrastructure. A conflict of interest exists: Jan Christoph is registered as an inventor of patent application PCT/EP2015/077001.

All source code used in this article is available in an online repository under https://gitlab.com/janlebert/mechano-electric-synchronization.

The simulation and synchronization parameters for the four reconstructions shown in Figs. 4–7 are listed in Table I. Figure 13 depicts the integration scheme of the numerical simulation. The different steps for the generation of synthetic deformation data used for the reconstruction are shown in Fig. 14.

In experiments, we imaged ventricular tachycardia and fibrillation with approximately 80–160 volumes per second and fields of views ranging in the order of $50×50×50$ to $100×100×100$ voxels using high-resolution ultrasound.²⁸ With dominant frequencies of 5–10 Hz, which are typically seen during ventricular tachycardia or fibrillation, one rotation of a spiral wave rotor takes place in about 100–200 ms. With temporal resolutions of approximately 6–12 ms, such a rotation can be resolved with about 15–30 images. We anticipate that higher volume rates during ultrasound imaging will be feasible in the future.

Phase of the excitation was obtained using the Hilbert transformation. To avoid artifacts in the phase maps during the coupling of system 2 with system 1, the phase was derived only from every $τ$th time step, during which the coupling and synchronization were minimal. Phase singularities were computed as nonvanishing sums of the circular integral summing over the gradient of the phase around points $r$ in the phase maps,⁹