Characterization of multiple spiral wave dynamics as a stochastic predator-prey system

A perspective on systems containing many action potential waves that, individually, are prone to spiral wave breakup is proposed. The perspective is based on two quantities, “predator” and “prey,” which we define as the fraction of the system in the excited state and in the excitable but unexcited state, respectively. These quantities exhibited a number of properties in both simulations and fibrillating canine cardiac tissue that were found to be consistent with a proposed theory that assumes the existence of regions we call “domains of influence,” each of which is associated with the activity of one action potential wave. The properties include (i) a propensity to rotate in phase space in the same sense as would be predicted by the standard Volterra-Lotka predator-prey equations, (ii) temporal behavior ranging from near periodic oscillation at a frequency close to the spiral wave rotation frequency (“type-1” behavior) to more complex oscillatory behavior whose power spectrum is composed of a range of frequencies both above and, especially, below the spiral wave rotation frequency (“type-2” behavior), and (iii) a strong positive correlation between the periods and amplitudes of the oscillations of these quantities. In particular, a rapid measure of the amplitude was found to scale consistently as the square root of the period in data taken from both simulations and optical mapping experiments. Global quantities such as predator and prey thus appear to be useful in the study of multiple spiral wave systems, facilitating the posing of new questions, which in turn may help to provide greater understanding of clinically important phenomena such as ventricular fibrillation.

Figure 1

(Color online) Definitions of the predator and prey states as applied to the time record of a single pixel taken from typical filtered optical mapping data. The pixel was considered to be in the predator or prey state when the optical signal, assumed to represent the membrane potential, was above or below fixed cutoff levels. Here the data were normalized and inverted so that 0.0 and 1.0 represented the minimum and maximum observed membrane potentials, respectively, for the given pixel. The predator and prey quantities $D\left(t\right)$ and $Y\left(t\right)$, at any given time $t$, were then defined to be the fraction of cells in the predator and prey states at time $t$.

Figure 2

(Color online) (a) Predator and prey quantities [$D\left(t\right)$ and $Y\left(t\right)$] vs time during type-1 behavior in the 3V-SIM system. The spiral wave period $\left(101\mathrm{ms}\right)$ is shown for reference. (b) Snapshots of the membrane potential (resting $\text{state}=0$, activated state $\approx 1$) in space taken at various times during the time interval shown in (a). The first frame also shows the sense of rotation of the tips of the eight spiral waves present. All three frames in the first row illustrate these rotations by marking the locations of the eight spiral wave tips (identified by the numbers 1–8). Frames in the second row, which are taken $100\mathrm{ms}$, or almost exactly one spiral wave period $(=101\mathrm{ms})$, later than those in the first row, show the continued, nearly identical pattern of rotation of these eight spiral waves.

Figure 3

(Color online) (a) Predator and prey quantities [$D\left(t\right)$ and $Y\left(t\right)$] vs time during type-2 behavior in the 3V-SIM system. The time interval used is equal in duration to that used in Fig. 2a, to allow direct comparison. The spiral wave period $\left(101\mathrm{ms}\right)$ is shown for reference. (b) Snapshots of the membrane potential (resting $\text{state}=0$, activated state $\approx$ 1) in space taken at various times during the time interval shown in (a). Frames in the second row are taken $100\mathrm{ms}$, or almost exactly one spiral wave period $(=101\mathrm{ms})$, later than those in the first row.

Figure 4

(Color online) Typical trajectories in predator-prey phase space during type-1 (purple) and type-2 (orange) behavior. The arrows show the dominant sense of rotation for the trajectories shown. The predator and prey data shown in Figs. 2a, 3a are short segments of the data used to produce the type-1 and type-2 trajectories, respectively, shown here.

Figure 5

(Color online) The triangularly shaped accessible region of predator-prey phase space may be divided up into two regions that are useful in classifying the types of predator-prey behavior observed.

Figure 6

(Color online) (a) Predator and prey histories for one of the longer 3V-SIM runs showing transitions back and forth between type-1 and type-2 behavior and (b) power spectrum of the predator quantity vs time for the same run. System size: $6\mathrm{cm}\times 6\mathrm{cm}$. ${\tau }_w^{-}=41\mathrm{ms}$.

Figure 7

(Color online) Mean time to extinction (in ms) as a function of current location in predator-prey phase space. The data used for the calculation were taken from the 48 3V-SIM simulations with identical parameters and randomized starts for the initial spiral wave, as described in the Methods section (Sec. 2). Regions I (smaller triangle) and II (larger triangle) from Fig. 5 are shown for reference.

Figure 8

(Color online) Fraction of time spent in various regions of predator-prey phase space (i.e., phase space density) plotted as a function of the predator and prey variables, calculated using the same data used for Fig. 7. The phase space density is normalized to the maximum phase space density.

Figure 9

(Color online) Behavior of the predator and prey quantities for APD restitution curves with different slopes. (a) The APD-DI restitution curves for 3V-SIM dynamics for five different values of the parameter ${\tau }_w^{-}$. A line of slope 1 is shown for reference. The remaining panels show the predator and prey histories and predator power spectrum vs time for three of these restitution curves: (b) ${\tau }_w^{-}=25\mathrm{ms}$, (c) ${\tau }_w^{-}=41\mathrm{ms}$, and (d) ${\tau }_w^{-}=120\mathrm{ms}$, from 15-by-$15\mathrm{cm}$ 3V-SIM simulations.

Figure 10

(Color online) Predator (red) and prey (green) histories and predator power spectrum vs time for (a) a shorter FMG run with system size $13.5\mathrm{cm}\times 13.5\mathrm{cm}$ and (b) a longer FMG run with system size $20.5\mathrm{cm}\times 20.5\mathrm{cm}$.

Figure 11

(Color online) (a) Definitions of instantaneous period and amplitude illustrated using a portion of experiment No. 1 data (see below). Asterisks mark the calculated local extrema of the data. (b) Optical mapping snapshots of the membrane potential $(\text{red}=\text{activated})$ in a right ventricular endocardial preparation taken at the two times marked as “I” and “II” in panel (a). (c) Instantaneous periods $T$ relative to the spiral wave period plotted vs the corresponding amplitudes $A$ of the predator quantity for two simulations (simulation No. 1, ${\tau }_w^{-}=25\mathrm{ms}$; simulation No. 2, ${\tau }_w^{-}=40\mathrm{ms}$) and two experiments (experiment No. 1, right ventricular epicardial preparation perfused with low-calcium $\left(0.125\mathrm{mM}\right)$ Tyrode solution; experiment No. 2, right ventricular endocardial preparation with blebbistatin $20\mu \mathrm{M}$ added to normal Tyrode solution to reduce mechanical motion. Slopes of the least-squares fit of $\log \left(T\right)$ vs $\log \left(A\right)$: simulation No. 1, 0.538; simulation No. 2, 0.538; experiment No. 1, 0.474; experiment No. 2, 0.509. The bottom three data sets were offset downwards by 1, 2, and 3 orders of magnitude from the simulation No. 1 data set, for clarity. (d) Distribution of data vs the instantaneous period (relative to the spiral wave period; logarithmic binning) for the same simulations and experiments as in (a). (e) Instantaneous period $T$ vs instantaneous amplitude $A$ of the predator quantity multiplied by the square root of the fraction of total area used to calculate the predator quantity. Three data sets are shown, all from the same simulation $({\tau }_w^{-}=25\mathrm{ms})$, with the predator quantity calculated from the entire simulation region (green circles), 1/9th of the region (black “+’s”), and 1/36th of the region (blue triangles). The log-log fits least-squares slopes were found to be 0.538, 0.539, and 0.535, respectively. In both (c) and (e), the solid lines depict the least-squares fit; the dashed lines mark unit standard deviation departures from each fit.

Figure 12

(Color online) An example of how the domains of influence might be drawn for a particular pattern of action potential propagation. The background color image is a snapshot of a large $(15\mathrm{cm}\times 15\mathrm{cm})$ 3V-SIM simulation $({\tau }_w^{-}=41\mathrm{ms})$ taken during type-2 behavior. $\text{Red}=\text{predator}$ (action potential) regions; $\text{green}=\text{prey}$ regions; $\text{blue}=\text{neither}$ predator nor prey. Black and white dots mark spiral wave tips rotating in the clockwise and counterclockwise directions, respectively. The black lines mark the boundaries of the domains of influence. Domains located near the edges of the system continue through the periodic boundaries to the opposite side of the system. Note that each domain contains exactly one spiral wave, except the one marked with the “X” which contained two spiral waves until they annihilated one another shortly before the time of this snapshot. The boundaries drawn are not definitive, but rather represent the authors’ best guess of an average of where the adjacent spiral waves had collided, will collide, and/or would have collided had not breakup or some other disrupting event occurred.

Figure 13

Changes in the areas of the predator and prey regions during a time interval $\Delta t$ occur on the boundaries of these regions. An area $F\Delta t$ is added to the predator region and subtracted from the prey region as the predator region (alias, the action potential) advances into the prey region (containing excitable tissue). Predator area $B\Delta t$ is lost on the trailing edge of the action potential, while prey area $R\Delta t$ is regained along its boundary with refractory tissue as the latter recovers. Here $F$, $B$, and $R$ are the rates at which the corresponding areas are gained or lost per unit time.

Figure 14

Illustration of the implications of wave backs and the leading edges of the prey region following and mimicking the propagation of wave fronts. The location of the wave front at time $t$ (represented by the thick line in all three panels) is the location of the wave back at the later time $t+\mathrm{APD}$ [panel (b)], which is also the location of the leading edge of the prey region at time $t+\mathrm{REC}$ [panel (c)]. Thus, the area swept out by the wave front between times $t$ and $t+\Delta t$ [$F\Delta t$ in panel (a)] is the same area swept out by the wave back between times $t+\mathrm{APD}$ and $t+\mathrm{APD}+\Delta t$ [$B\Delta t$ in panel (b)] which is the same as that swept out by the leading edge of the prey region between times $t+\mathrm{REC}$ and $t+\mathrm{REC}+\Delta t$ [panel (c)].

Figure 15

(Color online) (a) Qualitative illustration of the phase shift in time of the oscillatory part of the prey variable $Y\left(t\right)$ compared to that of the predator variable $D\left(t\right)$. (b) Phase space trajectory of $\mathbf{(}D\left(t\right),Y\left(t\right)\mathbf{)}$ when the peak in $Y\left(t\right)$ lags behind the minimum in $D\left(t\right)$.

Figure 16

(Color online) Data from a simple domain-based model. (a) Instantaneous periods $T$ relative to the spiral wave period plotted vs the corresponding amplitudes $A$ of the predator quantity for $N_d=30$. Calculated least-squares slope: 0.515. (b) Distribution of data vs the instantaneous period (relative to the spiral wave period, logarithmic binning) for the same model run as in (a). (c) Instantaneous period $T$ vs instantaneous amplitude $A$ of the predator quantity multiplied by the square root of $N_d$ for three different values of $N_d$: $N_d=10$ (blue dots), 40 (black “+’s”), and 360 (green circles). The log-log fit least-squares slopes were found to be 0.534, 0.531, and 0.523, respectively. (d) Typical behavior of the global predator quantity as a function of time $t$ taken from the $N_d=30$ run. In both (a) and (c), the solid lines depict the least-squares fit; the dashed lines mark unit standard deviation departures from each fit.