Data-driven phase-isostable reduction for optimal nonfeedback stabilization of cardiac alternans

Phase-isostable reduction is an emerging model reduction strategy that can be used to accurately replicate nonlinear behaviors in systems for which standard phase reduction techniques fail. In this work, we derive relationships between the cycle-to-cycle variance of the reduced isostable coordinates for systems subject to both additive white noise and periodic stimulation. Using this information, we propose a data-driven technique for inferring nonlinear terms of the phase-isostable coordinate reduction framework. We apply the proposed model inference strategy to the biologically motivated problem of eliminating cardiac alternans, an arrhythmia that is widely considered to be a precursor to more deadly cardiac arrhythmias. Using this strategy, by simply measuring a series of action potential durations in response to periodic stimulation, we are able to identify energy-optimal, nonfeedback control inputs to stabilize a period-1, alternans-free solution.

Figure 1

Steady state pacing of the Noble model from Appendix pp1. Panel (A) illustrates representative period-1 behavior that occurs for slower pacing rates. Panel (B) shows representative period-2 behavior that emerges for faster pacing rates. This stable period-2 behavior is known as action potential duration alternans. On each beat, action potentials are taken to be the time between successive crossings of the $-70$ mV threshold, as denoted above with a dashed horizontal line.

Figure 2

The panels above show relevent terms from the reduction Eq. (6). The Noble model from Appendix pp1 is considered with a pacing rate of 290 ms. As described in the text, two action potentials are used to define the periodic orbit so that the resulting Floquet multipliers are positive. Panel (A) shows the voltage trace on the periodic orbit, with panels (B) and (C) giving the corresponding values of $I\left(\theta \right)$ and $C\left(\theta \right)$, respectively.

Figure 3

Numerical confirmation of the derived relationships between changing Floquet multipliers and the variance of the resulting isostable coordinates for simulations of Eqs. (16) and (23). Black dots in panel (A) show how the variance of the isostable coordinates changes when input $u\left(t\right)=\alpha sin(8\pi t/T)$ is applied. Predictions computed according to Eq. (25) are shown as open circles. Histograms of the measurements of the isostable coordinates obtained from simulations of Eqs. (16) and (23) are shown in Panels (C) and (D) providing a visual representation of the change in variance that results when sinusoidal input is applied. Note that overlapping regions of the histograms in panels (C) and (D) appear brown.

Figure 4

While the APD, defined in Eq. (28) is a more natural experimental definition, both $t_1^m$ and $t_2^m$ are random variables, which would complicate the analysis below. By contrast, the AAPD defined in Eq. (29) consists of only one random variable (since $t_0^m$ is deterministic). In practice, since $\mathrm{var}(t_1^m-t_0^m)$ is small, the difference between the variance of the APDs and AAPDs is negligable. The voltage threshold, $\mathrm{\Pi }$, defines a Poincaré section that is used to calculate APDs.

Figure 5

Illustration of the proposed strategy for estimating $C\left(\theta \right)$ of the reduced model Eq. (6) for a stable entrained orbit. The gray lines in each panel show the exact value of $C\left(\theta \right)$ calculated directly from the model equations using methods described in Ref. [18]. Black lines in each panel show the data-driven approximation of $C\left(\theta \right)/\mu =a_0/2+{\sum }_{k=0}^n\left[b_{k}sin(2\pi kt/T)+a_{n}cos(2\pi kt/T)\right]$ for the indicated value of $n$. The coefficients of the Fourier series expansion are determined using the strategy detailed in Sec. 3e1. Dashed lines show the first $n$ Fourier modes of the exact solution. The magnitudes of the resulting plots are appropriately normalized to provide visual comparisons with the exact solutions. While the proposed estimation strategy is unable to capture all of the details associated with the high frequency modes, the contributions of the lower frequency modes are accurately estimated.

Figure 6

Various periodic inputs $u\left(t\right)$ (in $\mu \mathrm{A}/{\mathrm{cm}}^2)$ are applied [panels (A), (C), and (E)] to the Noble model from Eq. (A1) using a pacing rate of 253 ms. When no input is applied, alternans emerges in steady state as seen in the histogram in panel (B). Occasionally, after a particularly large action potential duration, the ensuing action potential fails, leading to an APD of approximately 350 ms immediately afterward. The application of the sinusoidal inputs in panels (C) and (E) stabilize the unstable period-1 alternans-free orbit, thereby eliminating alternans as can be seen in corresponding panels (D) and (F), respectively.

Figure 7

Green lines in panels (A) and (C) represent the baseline stabilizing stimulus $u_1\left(t\right)$ (in $\mu \mathrm{A}/{\mathrm{cm}}^2)$ as described in the procedure from Sec. 3e2. Green bars in panels (B) and (D) show the respective histograms of the APDs. Black, red, and blue lines in panels (A) and (C) show stimuli $u_2\left(t\right)-u_1\left(t\right)$, highlighting the additional inputs added to $u_1\left(t\right)$. Histograms of corresponding color are shown in panels (B) and (D). Differences between the resulting variances of the measured action potentials are then used to infer information about the shape of $C\left(\theta \right)$.

Figure 8

Results from implementing the proposed strategy from Sec. 3e2 for estimating $C\left(\theta \right)$ from the reduced equations of the form Eq. (6) associated with an unstable periodic orbit. Each panel shows 100 independent estimates of $C\left(\theta \right)$ (colored lines) using the indicated number of APDs to estimate the required APD variances according to the Fourier series expansion $C\left(\theta \right)/\mu =a_0/2+{\sum }_{k=0}^n\left[b_{k}sin(2\pi kt/T)+a_{n}cos(2\pi kt/T)\right]$. The magnitudes of the resulting plots are appropriately normalized to provide visual comparisons with the exact solutions. The black lines in each panel show the exact value of $C\left(\theta \right)$ calculated numerically from the model equations using methods described in Ref. [18]. Taking more terms of the Fourier series expansion (i.e., with $n$ larger than 6) only yields small differences in the resulting values of $C\left(\theta \right)$ and these results are not shown.

Figure 9

Using the proposed strategy to infer the curve $C\left(\theta \right)$ with Fourier modes up to the indicated value of $n$ [and also using the exact value of $C\left(\theta \right)]$, the right panel shows the resulting optimal stimulus (in $\mu \mathrm{A}/{\mathrm{cm}}^2)$ that stabilizes alternans and yields a variance of 66 $\pm 1$ ms in the APDs. Corresponding bars in the left panel represent value of ${\int }_0^{T_{\mathrm{pace}}}{u}^2\left(t\right)dt$ for the indicated stimuli. Resulting energy consumption using the feedback (FB) control strategy from Eq. (39) taking $\gamma =-0.025$ is also shown, where the reported energy usage represents the average per pacing cycle. The input using the true value of $C\left(\theta \right)$ is optimal when using the nonfeedback strategy and the approximations that result when using the inferred curves provide better estimates of the optimal inputs as more terms of the Fourier series expansion are included. While the delayed feedback control uses less energy than the proposed nonfeedback strategy, practical drawbacks associated with the implementation of feedback control in experimental situations may outweigh the benefit from the energy savings. As such, the delayed feedback control strategy and the nonfeedback control strategy cannot be directly compared solely on the basis of overall energy usage.