Spiral-wave dynamics in excitable media: Insights from dynamic mode decomposition

We use two mathematical models for cardiac tissue to illustrate the application of DMD to the study of spiral waves in excitable media. The first is the set of two-variable coupled partial differential equations (PDEs) called the Barkley model [27]; and the second is the biophysically realistic O’Hara-Rudy (ORd) model [28]. The Barkley-model PDEs are:

$u$ and $v$ are the fast-excitation and slow-recovery variables, respectively, at the point $\bf{r}$ and time $t$; the time-scale separation between $u$ and $v$ is controlled by the value of $\epsilon$; the parameter $a$ sets the duration of excitation, and $\frac{b}{a}$ sets the threshold of excitation. $D$ is the diffusion constant of the medium (we use $D=1$). We solve Eqs. (1) and ( 2) by using the forward-Euler method for time marching and a five-point stencil for the Laplacian. The temporal and spatial resolutions are set to be $\Delta x=0.2$ (arbitrary units $au$) and $\Delta t=0.005$ $au$.

For our study with inexcitable obstacles in the medium, we use the O’Hara-Rudy (ORd) model [28] for cardiac myocytes. In a homogeneous medium, the ORd model uses the following PDE for the transmembrane potemtial $V_{m}(\bf{r},t)$:

here, $C_{m}$ is the membrane capacitance, $D$ the diffusion coefficient (for simplicity, chosen to be a scalar), and the total ionic current $I_{\rm{ion}}$ is a sum of the following ion-channel currents: the fast inward Na⁺ current $I_{Na}$; the transient outward K⁺ current $I_{to}$ ; the L-type Ca²⁺ current $I_{CaL}$; the Na⁺ current through the L-type Ca²⁺ channel $I_{CaNa}$ ; the $K^{+}$ current through the L-type Ca²⁺ channel $I_{CaK}$ ; the rapid delayed rectifier $K^{+}$ current $I_{Kr}$; the slow delayed rectifier $K^{+}$ current $I_{Ks}$; the inward rectifier $K^{+}$ current $I_{K1}$ ; the Na${}^{+}/$Ca²⁺ exchange current $I_{NaCa}$; the Na${}^{+}/$K⁺ ATPase current $I_{NaK}$; the Na⁺ background current $I_{Nab}$; the Ca²⁺ background current $I_{Cab}$; the K⁺ background current $I_{Kb}$; the sarcolemmal Ca²⁺ pump current $I_{pCa}$; for a full list of these currents and the equations that govern their evolution we refer the reader to Refs. [28, 29, 30], which also describe the finite-difference numerical methods that we use. In both the models we study, we restrict ourselves to two spatial dimensions and we employ no-flux boundary conditions.

Excitation waves in the Barkley model [Eqs. (1) and ( 2)] have small wavelengths, so they are vulnerable to wavebreaks, especially in the presence of heterogeneities in the medium. In contrast, the ORd model [Eq. (3)] yields waves with large wavelengths; these waves are more stable in heterogeneous media than their Barkley-model counterparts. As our objective is to investigate the detection of the phase singularity of a stable spiral wave, we use the ORd model for our study with a heterogeneous medium; the coupling between cells in the presence of inexcitable obstacles is modelled as in Refs. [31, 30].

The tracking of the tip of a spiral wave, as described in Ref. [13], is based on the idea that the normal velocity of the spiral wave at its tip is zero. We illustrate this for the Barkley model Eqs. 1-2, in which $du/dt=0$ at the spiral tip. We choose isopotential lines with value $u_{iso}=0.4-0.5$; we then track the intersection point of $u(x,y)$ with $u_{iso}$ at different times, i.e., we record the position $(x,y)$, at a given time, where.

The locus of these intersection points gives the spiral TT.

The data-driven DMD method employs a linear operator to model the spatiotemporal evolution of fields in a complex, typically nonlinear, system. If $\mathbf{x}_{n}$ represents the vector form of some spatial data (typically an image) at the $n^{th}$ time instant, then the best estimate of the linear operator $\mathcal{L}$ that translates $\mathbf{x}_{n}$ to its value $\mathbf{x}_{n+1}$, at the next time step, follows from the minimization problem

where $\mathinner{\!\left\lVert.\right\rVert}_{2}$ represents the $L^{2}$-norm, and ${m}$ is the total number of images (spatial data), each one of which is separated from its predecessor and successor by the sampling-time interval $\tau$. The eigenvalues and eigenvectors of $\mathcal{L}$ contain useful information about the evolution of the system and can be calculated efficiently by using singular value decomposition [32, 23]. In our analysis we choose ${m}$ such that it is greater than all the macroscopic time-scales (spiral-rotational periods) that are present in our system. We construct a matrix $X_{1}$, with column vectors of images at discrete times labelled $0,1,\ldots,(m-1)$, and a similar matrix $X_{2}$, with column vectors of images at discrete times labelled $1,2,\ldots,m$ as follows:

The operator that best fits Eq. (5) is

where $X_{1}^{\dagger}$ is the Moore-Penrose pseudoinverse [33, 34]. By using the DMD algorithm [23], we get the following spectral decomposition:

here, $\lambda_{i}$ and $\Phi_{i}$ denote the eigenvalue and eigenmode of $\mathcal{L}$, respectively. Depending on whether $|\lambda_{i}|$ $>1$, $|\lambda_{i}|$ $=1$, or $|\lambda_{i}|$ $<1$, the corresponding eigenmode $\Phi_{i}$ grows, remains constant, or decays, respectively, in time. In Section I of the Supplemental Material [35] we give the SVD-based method that we use to get a low-rank version of $\mathcal{L}$ and thence the dominant eigenmodes. We refer the reader to [23] for a detailed discussion of DMD.

We characterize spiral-wave dynamics here by spiral-TT patterns and the wave-rotation frequencies. We illustrate this for the Barkley model [Eqs. (1) and ( 2)]. To obtain the frequencies of a spiral wave, we record the time series of $u$, Fig. 1 (B), from a representative point, which is marked by a yellow star in the simulation domain [Fig. 1 (A)]. From the power spectrum $E(\omega)$ of this time series we obtain the most important frequencies, like the frequency $\omega_{D}$ of the highest peak in this spectrum. We then use the IIM and DMD methods to obtain the TTs. In the illustrative plots of Fig. 2 we show four different types of TTs, in each one of the panels (A)-(D); the sub-figures in the left columns show the TTs obtained from the IIM; those in the middle columns show, via pseudocolor plots of the modulus of a DMD eigenmode $\Phi$, with eigen value $|\lambda|=1$, that contains an imprint of the TT; the right columns show the power spectra $E(\omega)$. The circular TT in Fig. 2 (A) arises from periodic motion, with a fundamental frequency $\omega_{D}=0.05$, which appears clearly in the power spectrum along with its harmonics. The TTs in Figs. 2 (B), (C), and (D) exhibit meandering patterns with petals on a circular trajectory, rosettes, and a rectangular trajectory with petals, respectively; these patterns arise from aperiodic motions with more than one fundamental frequency; e.g., the peaks in the $E(\omega)$ in Fig. 2 (B) can be labelled as $n_{1}\omega_{D}+n_{2}\omega_{2}$, with $n_{1}$ and $n_{2}$ integers (positive or negative), $\omega_{D}=0.253333$ and $\omega_{2}=0.245926$, and the ratio $\omega_{2}/\omega_{D}\simeq 0.9707$ an irrational number. We show other TT patterns, for different parameter sets, in Fig. S2 in the Supplementary Material [35]. From Fig. 2 we conclude that both IIM and DMD methods can detect complicated TT patterns equally well in a homogeneous medium.

We show, in Fig. 3, how the TT patterns vary with the non-dimensionalised sampling interval $\tau=\mathcal{T}\times\Delta t\times\omega_{D}$ (Table 1) in both the IIM and DMD methods. For specificity, we use the rosette pattern in Fig. 3 (C). The TT pattern is roughly conserved as we increase $\tau$; however, the IIM yields noisy TTs for $\tau\geq 0.015$, whereas the DMD method yields a TT pattern that is less noisy than its IIM counterpart.

We now demonstrate that DMD can be used to determine the positions of phase singularities even when there are multiple spiral waves in the medium. For the Barkley model [Eqs. (1) and ( 2)], we illustrate this in Fig. 4 (A), which shows a pseudocolor plot of $u$ in a multiple-spiral state; and Fig. 4 (B) depicts the modulus of a DMD eigenmode $\Phi$ that exhibits clearly the locations of the phase singularities that are present in this pseudocolor plot. In Fig. 4 (B) we also see the boundaries between the domains with spiral waves.

We use the O’Hara-Rudy model [28] of cardiac excitation waves for our study of TTs in a heterogeneous medium. Figure 5 shows how DMD and IIMs track TT in the medium with two different percentages $P_{o}$ of obstacles. We find that, although both methods can locate the region where the TT is confined, the TT plot from the IIM is associated with randomly scattered points in the background, which are suppressed in the TT pattern extracted by the DMD eigenmode. Moreover, we check how these two methods perform in the presence of noise in the signal. Such noise can arise in the data-collection processes in real experiments. Figure 6 shows TTs for three different values of signal-to-noise (SNR). It shows that IIM is more sensitive to noise and it fails to track the TT for SNR $<22$, whereas DMD can still capture the TT pattern. DMD can produce the TT pattern upto SNR$\simeq 16$. In summary, our results demonstrate that, with external noise, DMDTT is a more robust and versatile method for tracking TTs than IIM.

The DMD technique can also be used to reconstruct (approximately), and hence predict, the spatiotemporal evolution of the spiral waves in the models we consider. For this reconstruction we require the eigenmodes $\Phi_{i}$, their eigenvalues $\lambda_{i}$, and the amplitudes $b_{i}$ associated with every mode $\Phi_{i}$ (see, e.g., Ref. [25]). We illustrate this for the Barkley model [Eqs. (1) and (2)], where we predict the spiral-wave dynamics as follows; at any instant of time $t$, the predicted solution $u^{P}{(\bf{r},t)}$ is

where the coefficient $b_{i}=(\Phi^{{\dagger}}\Phi)^{-1}\Phi^{{\dagger}}({\bf{r}})x_{0}$, the ${\dagger}$ denotes Hermitian adjoint, $m$ is the number of columns in $X_{1}$, and $x_{0}$ is the first column in $X_{1}$ (Eq.( 6)). We show in Figs. 7(A), (B), and (C) pseudocolor plots of $u({\bf{r}},t)$ [from Eqs. (1) and (2)], the DMD prediction $u^{P}({\bf{r}},t)$ [from Eq. (9)], and the error $u({\bf{r}})-u^{P}({\bf{r}},t)$, respectively. From the plots in Figs. 7(A), (B), and (C) we conclude that the DMD-based spiral-wave reconstruction works well here because the error $u({\bf{r}},t)-u^{P}({\bf{r}},t)\simeq 10^{-2}$. We suggest that such a DMD-based prediction can be used mutatis mutandis with potentiometric voltage data for spiral waves from ex-vivo and in-vitro experiments.