Estimating the probabilities of rare arrhythmic events in multiscale computational models of cardiac cells and tissue

Fig 1. Multiscale cell and tissue model schematics. (a) Diagram of intracellular Ca2+ compartments and transport. At each release site at coordinates (i,j,k), Ca2+ is released via RyRs from the JSR into the dyadic subspace (SS) and diffuses into a submembrane (SM) compartment. Ca2+ diffuses between SM compartments of adjacent release sites in the 3D lattice, as depicted in panel (b). Ca2+ can also diffuse from the SM into a single cell-averaged cytosolic compartment and can be transported by the SR Ca2+ ATPase (SERCA pump) from the cytosolic compartment into a single network SR (NSR), which refills the JSR. (c) Illustration of the different spatial scales incorporated in the tissue model including (from left to right) the 48 stochastic RyRs in each of the 25,000 release sites of a single cell. Hundreds of cells are coupled via gap junction currents to form the fiber model.

Fig 2. Ventricular myocyte model properties. (a) Normalized peak LCC (triangles) and RyR (circles) flux when the membrane was stepped to test potentials between -20 and 40 mV. (b) Excitation-contraction coupling gain, defined as the ratio of the normalized peak RyR to LCC Ca2+ flux at each test potential. (c) Representative action potentials and (d) [Ca2+]i transients under control conditions (solid black), with 50% IK1 inhibition (dashed black), with β-adrenergic stimulation (solid red) conditions, and with β-adrenergic stimulation and 50% IK1 inhibition (dashed red). (e) Cell-wide SR Ca2+ leak rate via RyRs at varying SR Ca2+ loads in the baseline model. The dotted line indicates the lowest SR Ca2+ load tested that exhibited spontaneous Ca2+ waves. (f) Example linescan plot showing a Ca2+ wave under β-adrenergic stimulation.

Fig 3. DADs induced by Ca2+ overload during 1 Hz pacing under β-adrenergic stimulation in the myocyte model. (a) Sub-threshold DADs (red arrows) and triggered APs (blue arrows) in between paced APs. Inset illustrates how the first DAD occurred just prior to the next stimulus. (b) Spontaneous [Ca2+]i transients caused by underlying Ca2+ waves. Dotted lines indicate diastolic [Ca2+]i. (c) Average [Ca2+]JSR. Dotted lines indicate threshold for Ca2+ overload that resulted in spontaneous Ca2+ waves. This threshold decreased over the course of the protocol due to increased RyR sensitivity, increased [Na+]i, and elevated diastolic [Ca2+]i.

Fig 4. Elevating SR Ca2+ load accelerated Ca2+ wave formation and increased DAD amplitude. (a) DADs resulting from initializing SR Ca2+ load to five different values as shown in (b). (c) Volume renderings of Ca2+ in the simulations at three time points illustrating the greater spontaneous Ca2+ wave activity at higher SR Ca2+ loads. The number of Ca2+ wave nucleation sites (Nnuc) is tabulated for each simulation. All simulations initialized with [Na+]i = 10 mM.

Fig 5. Variability of DAD timing and amplitude in five independent model realizations with identical initial conditions. (a) Variability in DAD timing and amplitude. Inset shows the area under the curve above resting potential of each DAD. Strong correlations were observed between Vmax and the maxima of both the spontaneous [Ca2+]i transient (b) and the RyR Ca2+ release flux (c). (d) Volume renderings of Ca2+ wave dynamics at three time points in each simulation. Approximate locations of the nucleation sites are circled in white. Nnuc was higher in simulations i and ii which had greater DAD amplitude. All simulations initialized with [Na+]i = 10 mM.

Fig 6. Roles of cytosolic and SR Ca2+ and IK1 density in DAD distribution. DAD distributions when initial (a) [Ca2+]SR and (b) [Ca2+]i were varied in the presence of β-adrenergic stimulation. (c) DAD distributions when varying initial [Ca2+]i after reducing IK1 density by 50%. Shaded regions indicate the middle quartiles of the membrane potential for sub-threshold events out of 96 simulations. Dark lines indicate median values. Statistics of Vmax (d) of sub-threshold events, phase plane plot of V and NCX current (e), and statistics of delay until DAD peak (f) are shown when SR [Ca2+] was varied (blue), [Ca2+]i was varied (red), and [Ca2+]i was varied with 50% reduction in IK1 density (magenta) as a function of initial total cell Ca2+. Error bars indicate standard deviation (SD) of the estimates. Panel (e) shows representative sub- and supra-threshold DADs for the baseline model and IK1 reduced by 50% with [Ca2+]i initially at 3.5 and 2 μM, respectively. Arrows indicate direction of time. (g) Dependence of triggered AP probability on initial total cell Ca2+. All simulations initialized with [Na+]i = 10 mM.

Fig 7. Probabilistic triggered activity caused by DADs in the fiber model. (a) Fiber of 192 cells paced at 0.5 Hz by applying a stimulus current for 2 ms to the first two cells in the fiber (black arrows). DADs occurred between paced beats (red arrows). (b) Ensemble simulations (i)-(iv) initialized to the state at time trestart in panel (a). Realization (i) corresponds to the simulation in panel (a). Enhanced view of (iv) is shown at right.

Fig 8. Roles of initial [Ca2+]i, IK1 density, and ggap in DAD variability in the fiber model. (a) Spatiotemporal profile of V (left) and Vmax profile (right) in a 480-cell fiber under baseline conditions. Similar simulations are shown with 50% IK1 reduction (b) and both 50% IK1 and 50% ggap (c). The range of Vmax values are indicated for the red traces, which correspond to the images on the left. To avoid boundary effects, the outer 24 cells on either end were initialized to normal SR loads. The inner cells were initialized to identical initial conditions to those of Fig 6. All simulations initialized with [Na+]i = 10 mM.

Fig 9. Filtering method for estimating Vmax from the spatiotemporal JRyR profile in a fiber. (a) Illustration of the filtering method for an example fiber simulation. A uniformly-weighted spatial smoothing filter of width W was applied to the spatiotemporal JRyR profile. The maximum values of the filtered profile were then normalized to obtain an estimate of the voltage profile (see Eq 6). (b) The filtering method was applied to three fiber simulations using the conditions indicated above each plot. The filtering method estimate (red) accurately reproduced Vmax from the original simulation (black). Parameters are listed in Table 1.

Fig 10. Method for estimation of rare extreme DAD probabilities. (a) Example illustrating how independent fiber realizations are generated by shuffling cell positions and applying the filtering method. The resulting Vmax profile of the original simulation (black) and Vfmax obtained after shuffling (red) are plotted on the right. Vpeak is defined as the maximum potential achieved along the length of the fiber. (b) Histograms of Vpeak from model simulations (black), the filtering-only method (blue), and filtering and shuffling method (red). See text for details. (c) Upper tail of the distributions from (b). (d) Predicted tail of Vpeak distributions plotted relative to μV for fibers with the baseline model (blue), with 50% IK1 reduction (red), and with both 50% IK1 reduction and 50% ggap (magenta).

Fig 11. Realizations of rare extreme DADs. (a) Plots of Vmax—μV in the realizations from Fig 10D containing the most extreme DADs out of 106 trials. Curves correspond to the baseline model (blue), 50% reduction of IK1 (red), and both 50% IK1 and 50% ggap (magenta). (b) Spatiotemporal JRyR profiles in each realization from panel (a). Brackets and dotted lines mark the location of the filter window centered on the extreme DAD.