Right ventricle wedges from canine were prepared according to the experimental protocols approved by the Institutional Animal Care and Use Committee of the Center for Animal Resources and Education at Cornell University. Fluorescence optical mapping recordings of the membrane potential were recorded at a spatial resolution of 600 × 600 μm² per pixel for a grid size of 7.7 × 7.7 cm² and a temporal resolution of 2 ms at physiological conditions. For details of the experimental setup information, we refer to the previous studies (Fenton et al., 2009; Luther et al., 2011; Gizzi et al., 2013; Gizzi et al., 2017).

We based our numerical simulations on the four-variable minimal model for ventricular action potentials (Bueno-Orovio et al., 2008), solved on two-dimensional anisotropic heterogeneous spatial domains according to the fine-tuning performed by Fenton et al. (2013). The model includes a phenomenological description of main transmembrane ion currents and is properly generalized with a heterogeneous diffusion contribution to account for spatial effects. The model’s equations are ∂ t u = ∇ ⋅ D i j ∇ u − J f i + J s o + J s i , (3a) d t v = 1 − Θ u − θ v v ∞ − v / τ v − − Θ u − θ v v / τ v + , (3b) d t w = 1 − Θ u − θ w w ∞ − w / τ w − − Θ u − θ w w / τ w + , (3c) d t s = 1 + tanh k s u − u s / 2 − s / τ s , (3d) where u is the dimensionless cell membrane potential, v, w, and s are gating variables regulating ion current activation, and Θ(x) denotes the Heaviside step function. J _fi , J _so , J _si represent fast-inward, slow-outward, and slow-inward transmembrane currents: J f i = − Θ u − θ v u − θ v u u − u v / τ f i , (4a) J s o = 1 − Θ u − θ w u − u o / τ o + Θ u − θ w / τ s o , (4b) J s i = − Θ u − θ w w s / τ s i . (4c) Time constants and asymptotic values for gating variables depend on the membrane voltage: τ v − u = 1 − Θ u − θ v − τ v 1 − + Θ u − θ v − τ v 2 − , (5a) τ w + u = τ w 1 + + τ w 2 + − τ w 1 + tanh k w + u − u w + + 1 / 2 , (5b) τ w − u = τ w 1 − + τ w 2 − − τ w 1 − tanh k w − u − u w − + 1 / 2 , (5c) τ s o u = τ s o 1 + τ s o 2 − τ s o 1 tanh k s o u − u s o + 1 / 2 , (5d) τ s u = 1 − Θ u − θ w τ s 1 + Θ u − θ w τ s 2 , (5e) τ o u = 1 − Θ u − θ o τ o 1 + Θ u − θ o τ o 2 , (5f) v ∞ = 1 − Θ u − θ v − , (5g) w ∞ = 1 − Θ u − θ o 1 − u / τ w ∞ + Θ u − θ o w ∞ * . (5h) D _ij is the two-dimensional diffusion tensor, defined as D i j ≡ σ i j S o C m = D 11 D 12 D 21 D 22 , (6) where σ _ij is the conductivity tensor, S _o represents the cell surface-to-volume ratio, and C _m is the membrane capacitance. The tensor elements are defined in the two-dimensional Cartesian domain as D 11 = D ∥ ⁡ cos 2 α x , y + D ⊥ ⁡ sin 2 α x , y , (7a) D 22 = D ∥ ⁡ sin 2 α x , y + D ⊥ ⁡ cos 2 α x , y , (7b) D 12 = D 21 = D ∥ − D ⊥ cos α x , y sin α x , y . (7c) Here, the function α(x, y) represents the local fiber orientation and D _∥ and D _⊥ denote diffusivities along the directions parallel and perpendicular to the fibers. We used the same anisotropy settings as in previous computational studies on cardiac activation maps (Loppini et al., 2019). Model parameters are reported in Table 1. The parameter set is consistent with the one originally fine-tuned by Fenton et al. (2013), related to endocardial tissue at 37° Celsius. Specifically, parameter values are set to reproduce AP features, conduction velocity, and restitution curves as observed in canine cardiac tissues, the same considered in this study. Furthermore, we assume the anisotropy ratio as 1:3 in line with previous modeling analyses (Loppini et al., 2019). For a more detailed description of model parameters and for a comparison between the four-variable model results and experiments, we refer the reader to the two abovementioned studies.

As detailed in Section 3.4, the generalization to heterogeneous modeling by data assimilation is obtained by imposing a spatial variation of selected parameters, opportunely sorted around their reference values (i.e., Table 1) based on experimentally informed profiles. On this basis, the heterogeneous model is obtained by perturbing parameters of the homogeneous model so that global features in the evoked electrical activity are still correctly reproduced. We numerically integrated the model with an explicit Euler scheme implemented in Fortran, discretizing the spatial operators to account for heterogeneous diffusion and phase-field boundary conditions. We solved the model in both 1D and 2D domains, assuming zero-flux boundary conditions. Space and time discretization is Δx = 0.025 cm, Δt = 0.01 ms. Stability and conduction velocity convergence was verified upon mesh refinement and testing higher-order discretization schemes in time (second- and fourth-order Runge–Kutta), achieving non-significant variations in the computed results.

Simultaneous recordings of the epicardium and endocardium in RV canine preparations were observed (Gizzi et al., 2013). Physiological alternans-free wave conduction and alternans states could be observed depending on the pacing site position and pacing frequency. At lower pacing frequencies, no alternans is observed, as the APD is sufficiently short to prevent the interaction of subsequent waves. Figure 1A shows such an example observed at the entrainment frequency f _p = 3.2 Hz on the epicardium. The phase and amplitude of the pixel-wise FFI-analyzed recordings show a continuously evolving phase and amplitude in the entire tissue. No spatially correlated phase or amplitude is observed at f _1/2 = f _p/2 = 1.6 Hz that would indicate alternans. Also, closer inspection of the local signaling does not indicate alternans (Figure 1A, right panels). The normalized APs at P1 (pink square, top) and P2 (cyan square, bottom) show no alteration in peak height or APD, as confirmed by the respective amplitude in the Fourier space. Only a single peak at the entrainment frequency f _p is observed. The lower amplitude peak at 2f _p shows the second mode and does not carry relevant information. Contrarily, alternans is observed when the epicardium is entrained with a higher frequency. Figure 1B shows the same analysis and local signaling recordings at f _p = 8.0 Hz. In this case, the pacing frequency is sufficiently high so that subsequent waves influence each other. The phase and amplitude information at f _1/2 = 4 Hz shows a typical pattern of discordant 2:2 alternans. That means that two subsequent waves lead to two different APDs in time. In space, the APD of each wave can transiently switch between the two APDs that are spatially confined by nodal lines. The latter can be identified by the spatial phase jump of about π, and the amplitude valley. The normalized APs at P1 (cyan) and P2 (pink) show APD alternation between shorter and longer APDs. Every second AP is shown in either black or red to facilitate visualization. For those time series, a second peak at f _1/2 = 4 Hz is visible in the amplitude spectrum, since a second underlying frequency is present that correlates with two times the wave period (2T = f _1/2).

Although 2:2 alternans is the most commonly observed AP rhythm, other higher-order AP rhythms exist (see, e.g., Figure 6 in Gizzi et al. (2013)). The epicardium that is shown in Figure 1 was additionally paced from the bottom (base) of the heart at f _p = 8.5 Hz, which led to a spatially mixed mode of 2:2 and 4:4 alternans (Figure 2A). While the 2:2 AP rhythm shows a single amplitude peak at f _1/2 = 4.25 Hz, two additional amplitude peaks are observed in the Fourier space for the 4:4 AP rhythm: one very close to f _1/2 and one at f _3/4. As the two peaks at around f _1/2 are very close to each other but implicate different information, they are from here on defined as f 1 / 2 1 and f 1 / 2 2 . The well-defined spatial distribution of the two different AP rhythms is visualized at the phase and amplitude reconstructions in Figure 2B. Here, f 1 / 2 2 and f _3/4 show additional spatial information to the frequencies at f _p and f 1 / 2 1 in the Fourier space. The 2:2 (left side, RV anterior) and 4:4 (right side, RV anterior) alternans regimes are spatially separated at f _3/4. The right side shows synchronized phase and correspondingly elevated amplitude signaling that is present on the left side. At f 1 / 2 1 and f 1 / 2 2 , more complex phase dynamics are observed in addition to the nodal lines that separate different entrained alternating regions in the tissue. Figure 2C shows the corresponding phase wave directions. Numbers 1 and 2 indicate the temporal order of occurring phase waves. In this framework, nodal lines and nodal areas carry more information compared to the classically analyzed temporal differences in APD (Pastore et al., 1999). At the phase of f 1 / 2 1 (top, right side, Figure 2B), the phase wave accelerates at the nodal line from 1 to 2 (Figure 2C), while contrarily at f 1 / 2 2 , the phase waves propagate in the opposite direction and jump by 2π (top, right side, Figure 2B). The latter indicate the typical phase waves, as observed in the Fourier space ( f 1 / 2 2 ) at 2:2 AP rhythms (see also Figure 1). For the sake of completeness, Figure 2 shows additionally f _1/4 that corresponds to the duration of four AP rhythms. However, f _1/4 is a rather non-significant frequency and appears only because of the intrinsic noise of the system, as the peak at f _p describes already the main oscillation of the system, and the three other amplitude peaks at f 1 / 2 1 , f 1 / 2 2 , and f _3/4 describe the periodic temporal variations.

The frequency response observed at a single recorded pixel is useful to get an overview of the local alternans offset (in analogy to the well-known restitution curves). Figure 3A shows frequency maps with the normalized amplitudes depending on the entrainment f _p for top (base, left panels) and left (RV posterior, right panels) paced canine ventricles. The top and bottom panels show data recorded at the epicardium (EPI) and endocardium (ENDO). The main peaks (bright yellow peaks) indicate f _p. Above f _p appears a second peak from about 5 Hz that indicates alternans ( f 1 / 2 2 ) . Additionally, only for very high frequencies, here f _p = 9.2 Hz, a peak at f _3/4 is visible that indicates 4:4 alternans. Figure 3B shows a generalized scheme of the frequency maps for comparisons. The offset of 2:2 alternans differs between the EPI and the ENDO, while the peak f _3/4 is observed for all four frequency maps at the same entrainment frequency f _p = 9.2 Hz. At frequencies higher than 10.0 Hz, fibrillation is observed, as additionally shown in the right (f _p = 10.7 Hz, RV posterior) paced EPI and ENDO frequency maps. The comparison between f _p = 9.2 Hz and 10.7 Hz is shown in Figure 3C. Fibrillation shows an elevated baseline and a broad spectrum for frequencies above f _p.

While a critical frequency induces fibrillation, the complex spatiotemporal alternans patterns stabilize with the increasing entrainment frequency (Gizzi et al., 2013) (Figure 4A). Figure 4B shows selected snapshots of the EPI and ENDO from two different pacing sites. Initially, no alternans is observed at the EPI at a lower f _p ≃ 4 Hz, but the initiation of alternans at the ENDO can be seen (endocardium base paced, Figure 4). Interestingly, this occurs in larger speckled patches rather than in defined areas, which indicates the alternans-offset difference among individual cells. This speckled-like early fine-scale initiation of alternans was suggested previously by Jia et al. (2010). With increasing f _p, those patterns synchronize spatially and lead to distinct phase areas that are spatially separated by nodal lines, as best visible at the EPI. Although the ENDO shows comparable stabilization of alternans in the phase, the amplitude shows more spatial variations. This is most likely caused by the influence of the Purkinje fibers that are confined to the subendocardial layer and believed to be responsible for the initiation of ventricular fibrillation (Fox et al., 2002; Muñoz et al., 2018).

As the differences of the electrophysiological properties of individual cells also lead to differences in the alternans-offset and restitution characteristics, it is useful to take pixel-based differences into account when modeling alternans dynamics in silico. The advantage of the Fourier analysis of heart tissues is that the optical ultrastructure can be revealed in the low-frequency regime, as shown in Figure 5. Especially, the amplitude information at f = 0.5 Hz is a stable indicator for morphologically restricted differences that are independent of the pacing location (Figure 5A) and pacing frequency (Figure 5B). Here, we utilize the low-frequency regime, as it is also an indirect measure of the signal height, i.e., the observed baseline of the AP rhythms. Therefore, we assume that the strength of the emitted signal depends on the local tissue properties and thus relates to the heart ultrastructure.

In order to validate this hypothesis, we propose a novel data assimilation approach using the ultrastructure observed at f = 0.5 Hz assuming the influence in the diffusive term, Eq. 2.3 (D _∥, D _⊥), and the time constants that shape the AP, Eq. 5 ( τ w + , τ _so , and τ _si ). The tissue heterogeneities are applied on a pixel-based scale, so that the local variations of conduction velocity and action potential are accounted through specific parameters having a strong impact on the resulting APD. In detail, a mask encoding the actual tissue of the experimental samples was extracted by analyzing the signal-to-noise ratio, and the spatial map of the Fourier amplitude spectrum at f = 0.5 Hz was computed. The resulting ultrastructure profile was smoothed by using a Gaussian kernel, restricted to a radius of 6 pixels with a variance equal to 5. The heterogeneity field was evaluated by normalizing the ultrastructure mask as H x , y = δ r m a x | r | + 1 , (8) where r = s − s ̄ , s = |F _x,y (0.5 Hz)|, and δ is the parameter denoting the desired maximal variation. We set this value in the range [0,1] by assuming heterogeneity variations of the local parameters up to 100%. Eventually, the heterogeneity map is combined with the tissue map to include also information on tissue boundaries (i.e., phase-field). We refer to this heterogeneity field as an H ₁(x, y) map. Furthermore, we considered a second heterogeneity field that enhances the data assimilation procedure from low-amplitude areas in the Fourier spatial map. We evaluated the reciprocal of the H ₁ (x, y) map and normalized the result according to Eq. 8, still constraining the parameter δ ∈ [0, 1]. This second spatial scalar field is denoted as the H ₂(x, y) map with average 1 and variation δ. It is worth noting that the proposed procedure represents a generalization of the phase-field method (Fenton et al., 2005) that permits both the inclusion of tissue boundary and specific modulations of model parameters. The resulting heterogeneity maps are used as a multiplication factor on selected model parameters to achieve a constitutive heterogeneity shaped on a tissue ultrastructure. In our analyses, we focused on heterogeneities related to diffusivity and APD parameters (D _∥, D _⊥, τ w + , τ _so , τ _si ), by using the following constitutive law: p x , y = p ̄ H i x , y . (9) Here, p (x, y) denotes a spatial dependent parameter, p ̄ is the original model parameter, and H _i (x, y) is the heterogeneity field (with i = 1, 2). A visual representation of the adopted data assimilation technique is thoroughly provided in the next sections. Naming P ₁ as the set of diffusivity parameters and P ₂ as that of APD-regulating time constants, we investigated all possible combinations of H ₁ and H ₂ maps on P ₁ and P ₂, i.e., a specific heterogeneity profile was applied to diffusivities and APD-related time constants. It is worth noting that, with this setting, we are assuming a correlation between diffusivity and APD, with such parameters non-linearly coupled through complex electronic effects in cardiac tissue.

We performed an extensive in silico study on both one-dimensional cables and two-dimensional tissues to test the heterogeneity effects on alternans onset and severity. In particular, we computed H ₁ (x, y) and H ₂ (x, y) maps from a selected experimental tissue to shape model parameters’ heterogeneity in 1D and 2D domains, investigating all possible combinations of the heterogeneity fields on diffusivity and APD-regulating time constants, at δ = 0.25, 0.5, 1. For 1D simulations, we extracted one-dimensional cuts of H ₁ (x, y) and H ₂ (x, y) maps along the experimental propagating wavefronts (not shown). This preliminary set of numerical simulations was used as a first benchmark of the data assimilation procedure. We observed that the heterogeneous model is able to 1) recover the expected average CV and AP features and 2) emphasize alternans onset and severity, also inducing conduction block phenomena not observed in the homogeneous case. We then tested data assimilation within 2D computational domains observing notable differences with respect to the homogeneous case. In the following, we show two representative examples comparing the overall results for the same ventricle stimulated with a pacing-down protocol both at the ventricle base and in the right anterior ventricular region. The pacing-down protocol consists in stimulating the tissue starting from a low frequency and progressively reducing the pacing frequency. In particular, at each frequency, we delivered a stimulation train of 10 beats to ensure the tissue reached a stationary regime. This protocol reproduces the experimental one, and in our analyses, we computed alternans patterns on the last two beats to avoid transient effects.