A computational framework for testing arrhythmia marker sensitivities to model parameters in functionally calibrated populations of atrial cells Open Access

Cardiac tissue has long been recognized to behave as a highly complex dynamical non-linear system. The non-linear nature of the interactions between cellular and tissue components gives rise to characteristic pathological behaviors observed in the heart, such as cardiac alternans, spiral waves, and spatiotemporal chaos, ultimately leading to arrhythmia and fibrillation. Atrial and ventricular arrhythmias, in particular, are linked to multiple mechanisms interacting in a non-linear manner on a number of temporal and spatial scales, ranging from beat-to-beat to months and years, and from the subcellular to tissue levels.^1–4 On a subcellular and cellular scale, alterations of ion channel and pump distribution and function, intracellular ion dynamics, and metabolic pathways are the main causes of abnormal excitability.

One particular emergent behavior of interest for the study of cardiac arrhythmias is the adaptation of action potential (AP) to the pacing rate. This property is manifested at both cell and tissue levels and is often represented as the relationship between the AP duration (APD) and the previous diastolic interval or the pacing cycle length (PCL). Due to the so-called cardiac memory of the excitable membrane, rate adaptation and restitution properties are highly dependent on the initial conditions of the system. Depending on the pacing protocol chosen and on the particular electrophysiology (EP) of the cell (i.e., ion channel density and distribution), the APD restitution properties can either be stable, i.e., decreasing monotonically with faster pacing, or multistable.⁵ In the latter case, the cell can enter a period-doubling bifurcation, which is seen on the APD restitution curve as alternans, which are periodic alternations of long-short APs. The presence of APD alternans is generally seen as an indicator of cell instability, and has been shown to be often tightly linked to calcium transient alternans, i.e., beat-to-beat fluctuations of the calcium transient (CaT) amplitude from cycle to cycle.^4,6,7

Computational models constitute a powerful tool for studying the dynamics and underlying mechanisms of arrhythmia initiation and sustenance, by providing a quantitative framework to study the mechanisms underlying experimentally observed behaviors. Computer models have also been widely used to test new mechanistic hypotheses^8,9 and to predict putative clinical outcomes of pharmacological treatments, thus paving the way for in silico drug screening and therapy planning.¹⁰ 

Recent decades have shown numerous reports of experimental data on atrial EP, both in healthy and pathological conditions, and spanning from the single cell to organ. In parallel with the increasing amount of experimental data, a number of detailed and increasingly refined mathematical models have been developed for the atrial myocyte EP.^11–15 These models are typically constructed by fitting the model to averaged values of experimental measurements, with the aim of arriving at a single representative EP model. However, experimental findings have revealed a huge variability in the measured action potentials and ionic currents. These variations may result from different experimental conditions, i.e., “noise,” but may also be representative of inherent intra-subject or inter-subject variability in the atria electrophysiology. The sources of the observed variability may in part be attributed to differences in ion channel expression and dynamics, in particular the ionic currents I_CaL, I_to, I_Kur, I_K1, and I_NaK,¹⁶ which are among the most important currents in determining AP morphology.

In order to capture the variability observed in experimental data, an approach based on the so-called “population of models” has recently been proposed for the study of cellular electrophysiological mechanisms.¹⁷ A population of models generally refers to a set of models sharing the same ionic and molecular formulations, but with variable parameters to reflect observed variability in experimental measurements of the biomarkers. Conducting in silico experiments on populations of models allows the inclusion of inter- and intra-subject variability in the model system, and may uncover new emergent phenomena that are not observed in the traditionally single “averaged” model.^16,18,19

Sensitivity analysis (SA) is a technique that has been widely used in cardiac cell models to assess which parameters in the model exert the most influence on certain properties, and thus is particularly useful for studying dynamics of populations of models. In cardiac cell models, these properties are usually well-established EP markers, such as APD, AP amplitude, resting membrane potential, upstroke velocity, and afterdepolarizations.^16,20–22

A previous unpublished study using the Koivumäki model of the normal Sinus Rhythm (nSR) atrial cell²³ has shown evidence of unstable alternans behavior under dynamic pacing.²⁴ Simulations with the latest model version have shown the presence of APD alternans at unrealistic slow pacing rates compared to what has been observed experimentally. This bistable behavior is usually associated with chronic Atrial Fibrillation (cAF), as alternans in cAF-remodeled hearts are typically induced at much slower pacing rates as compared to nSR cells, and show impaired rate adaptation (i.e., flattened restitution).^25–27 This study aims to correct the physiologically unrealistic APD alternans behavior during dynamic pacing in this version of the Koivumäki human atrial cell model by fine-tuning specific model parameters. Additionally, we aim to elucidate the possible cellular mechanisms underlying APD alternans in both nSR and cAF by using a combined population of models and sensitivity analysis approach.

From the observation of the unstable behavior of the nSR model under dynamic pacing, we wanted to fine-tune the baseline model to mimic experimentally and clinically observed alternans behavior, while keeping the AP morphology observed in nSR. After fine-tuning the model to correct for the APD alternans behavior, we constructed and calibrated two populations of models representing both nSR and cAF, and performed sensitivity analyses on the populations to draw mechanistic conclusions regarding the influence of model parameters on pro-arrhythmia biomarkers.

We constructed a population of 1000 models using the Koivumäki atrial cell model (2016) in nSR²³ by varying 22 parameter values using a Gaussian distribution to randomly sample values around within ±30% of the baseline (σ = 0.15). The parameters that were varied and their baseline values are listed in Table I. The population of single cell models was implemented and solved using the single cell simulation environment BENCH (http://www.cardiosolv.com).

From the population, 14 single cell biomarkers were extracted, all based on properties of the AP morphology, and the APD alternans behavior. A list of the 14 biomarkers, and their definitions, is provided in Table II. All the chosen biomarkers are commonly used metrics to characterize single cell EP. The APD90, APD50, APD20, APtri, and APA are common metrics of AP morphology, which change characteristically in many pathological conditions, including cAF. The resting membrane potential (RMP) and dVdt_max are closely linked to the intracellular ion concentrations, and are important parameters of cell membrane repolarisation. The CaT amplitude, ΔCaT, was calculated as the difference between maximum and minimum (CaT₀) intracellular calcium concentration during one cycle. The time-to-peak of the CaT, CaT-tpeak, was determined from the onset of the CaT to the time of peak, and the time of decay, CaT-tdecay, was calculated as the interval from the peak to the point where CaT reached 10% of the CaT amplitude (Table III).

Biomarkers related to APD alternans were calculated from the APD (at 90%) restitution curve using a dynamic pacing protocol as follows: cells were first pre-paced at a PCL of 1000 ms for 30 s, and the PCL was decreased to 150 ms in steps of 10 ms. At each step, the cell was stimulated for 15 beats, and the last 5 beats were recorded for analysis. APD alternans were registered when the maximum and minimum APD values in the last 5 beats differed by more than 2%. Following the definition of the single cell biomarkers, chosen biomarkers were constrained within an acceptable range. This calibration step prevents model outliers that do not represent physiological behavior from being included in the sensitivity analysis. Choosing this calibration step is non-trivial, since the purpose of the population models is to capture variability, even beyond that observed in experimental results, but we wanted to exclude models that are obviously non-physiological. We calibrated the nSR population so that all models would fall within one standard deviation of the median values of APD90, APD50, APD20, APtri, APA, RMP, and dVdt_max. These biomarkers were chosen since they characterize AP morphology. Therefore, this type of calibration can be seen as “functional calibration,” in contrast to other calibration schemes employed.

In addition, we also excluded models with dynamic restitution curves exhibiting extreme maximum APD values. For the nSR population, the minimum and maximum APD accepted were 190 and 440 ms, respectively. These boundary values were taken from experimental reports of maximum and minimum APD90 values of nSR and cAF cells paced at 1000 ms PCL by Sanchez et al.¹⁶ 

From the nSR population generated in Sec. II A, we selected models with alternans threshold no larger than 250 ms (Alt-th ≤ 250). This particular value of alternans threshold was chosen for constraining the subpopulation since it is close to that reported by Narayan et al.²⁶ for patients in nSR. Imposing this constraint on the functionally calibrated nSR population resulted in subpopulation of models with no alternans or with alternans occurring only below 250 ms, whose mean parameter values were used to construct a new fine-tuned nSR baseline variant. This fine-tuned model was then used to represent the baseline nSR model in the rest of the study, according to the parameters listed in Table III.

We constructed two new populations of models, one representing nSR from the fine-tuned model created in Sec. II A and another representing cAF remodeling. cAF remodeling was simulated by the making the following changes on ionic components: 65% decrease in I_CaL, 18% decrease in I_Na, 68% increase in I_K1, 145% increase in I_Ks, 38% decrease in I_Kur, 50% decrease in I_KCa, and 62% decrease in I_to. The remodeling also accounted for a two-fold upregulation of RyR, 16% decrease in SERCA expression, two-fold increase in SERCA phosphorylation, and 10% cell dilation.

The two populations of models were generated and calibrated following the procedure described in Sec. II A. For the cAF population, the maximum and minimum accepted APD were 140 and 330 ms, respectively, as reported by Sanchéz et al.¹⁶ 

Then we performed sensitivity analysis on the functionally calibrated nSR and cAF populations in order to find linear correlations between responses (biomarkers) and multiple model predictors (model parameters). Sensitivity analysis was implemented based on the method proposed by Sobie et al.,²⁸ whereby chosen parameters from the model are simultaneously varied within a specified range to capture the “global” effect of a certain parameter on the model behavior. Sensitivity analyses were performed using the Multivariate Linear Regression (MVR) method.¹⁷ This method approximates non-linear dynamics to a linear model thus allowing us to identify key parameters that define the behaviors observed in the non-linear model. For each population, the varied model parameters were compiled into a matrix X of predictors, and the measured markers into a matrix Y of single cell biomarkers. X and Y matrices were centered and normalised prior to regression. The output of the MVR algorithm is a matrix B of regression coefficients calculated by (1)

where $Ŷ$ is an approximation of Y and ϵ is the regression error. The values b_ij in matrix B can be interpreted as surrogates of the sensitivity of marker y_j towards model parameter x_i. Larger values of b_ij indicate greater sensitivity, and positive b_ij means that an increase x_i is associated with an increase in y_j, while negative b_ij means that an increase in x_i results in a decrease in y_i. All analyses were implemented in Matlab2016a (https://mathworks.com).

In order to have a compact overview of the population sensitivities to the various parameters, we grouped the different cell biomarkers in AP-(APD90, APD50, APD20, APtri, APA, RMP, and dVdt_max), CaT-(ΔCaT, CaT₀, CaT-tpeak, and CaT-tdecay) and alternans- (Alt-th, Alt-range, and ΔAPD_max) related makers. For each group of biomarkers, we determined their overall sensitivity to the parameters by normalising the regression coefficients to the baseline and summing the contributions of the parameters to each marker in the group. We call these relative sensitivity values “global sensitivities.”

Figures 1(a) and 1(b) show the AP and CaT traces of the original nSR population at 1000 ms PCL. Figure 1(c) shows the dynamic restitution curves of the original nSR population before and after functional calibration, and Fig. 1(d) shows a few representative restitution curves for better visualization. The calibration step reduced the total population of 1000 models to 202 models.

As can be seen in Fig. 1(c), the original nSR population displayed a rather dynamic restitution behavior, with APD alternans occurring at all PCLs. The original nSR baseline model presents a single alternans region between 460 and 650 ms, with maximum ΔAPD of 21 ms. Figure 2 shows APD and calcium alternans in the nSR baseline model at a PCL 600 ms during dynamic pacing. The AP traces show a slight alternation in APD90 between long and short pulses of 18 ms (11% of APD90 of the long beat), with pronounced alternation of ΔCaT of ca 0.3 μM (60% of maximum ΔCaT).

Models from the calibrated population were then selected into a subpopulation of models showing APD alternans at PCLs equal or lower than 250 ms (Fig. 3). From this subpopulation, the mean parameter values were calculated in order to create a new baseline nSR model (nSR_tuned). This model did not exhibit APD alternans. However, the model did exhibit slight APD prolongation of 35 ms. The percent difference in the mean parameters relative to the original nSR baseline is shown in Fig. 4. The parameters that changed the most were iNaKmax and RyRtauinactss, and to a smaller extent, iCaPmax, RyRmaxss, and RyRtauactss. These parameters are related to alternans dynamics through the calcium cycling system. This results shows that a smaller change in these parameters is sufficient to alter the alternans behavior at the single cell level.

The simulated nSR and cAF populations are presented in Fig. 5. After functional calibration, the populations of 1000 models each were reduced to 213 models in the nSR and 357 models in the cAF population. The baseline, median, mean, and standard deviation values of the functionally calibrated nSR and cAF populations are shown in Table IV. The populations reproduced the AP morphologies characteristic of nSR and cAF cells. The nSR models presented a peak-and-dome shape, while the cAF models showed a more triangulated morphology, which is a known result of cAF ionic remodeling. Although the cAF baseline model did not present significant APD90 shortening as compared to the nSR baseline, which is a characteristic trace of cAF remodeling, the median APD90 of the functionally calibrated cAF population was shortened by 30 ms, and RMP hyperpolarized by 6 mV, as compared to nSR population. The CaT showed significant differences in nSR and cAF, with a 55% reduction of ΔCaT in cAF as compared to nSR, in agreement with previous studies.²⁹ Furthermore, while CaT-tpeak remained unchanged in cAF, CaT-tdecay was significantly longer in the cAF population as compared to the nSR population.

Figure 6 shows the distributions of APD90, APD50, APD20, APA, APtri, RMP, dVdt_max, ΔCaT, CaT₀, CaT-tpeak, CaT-tdecay, and Alt-th of the functionally calibrated nSR and cAF populations. The distributions of the biomarkers in the two populations presented in general different morphologies, with different mean values and standard deviations.

The sensitivities of APD90, APA, RMP, dVdt_max, ΔCaT, CaT_diast, Alt-th, and Alt-range are shown in Fig. 7. In our analyses, APD90 depends mostly on GCaL, GCab, GNa, GK1, and GKur, in line with previous reports.^{15,16,21,23,28} APD50 and APD20 showed different sensitivity trends in nSR and cAF. In nSR, APD50 is most dependent on GCaL, GKur, Gto, and GNa, while in cAF, APD50 depends on GKur, GNa, GCaL, and GNaL. APD20 in turn is most sensitive to Gto, GNa, GCab, and GK1 in nSR, and to Gto, GKur, GNa, and GCaL in cAF (not shown). Both APA and dVdt_max are most sensitive to GK1 and GCab, as well as GNa and INaKmax, with higher sensitivity in cAF. RMP is defined by GK1 and GCab in both nSR and cAF.

ΔCaT is primarily affected by GCaL, and to a smaller extent by Gto, RyRtauinactss, and kCa. CaT₀ is sensitive to Gab, INaKmax, and INaCamax. This is expected, since sarcolemmal calcium currents, in particular I_CaL, I_NaK, and I_NCX, are associated with intracellular calcium cycling. CaT-tpeak is also affected by GCaL, GCab, and GKur, RyRtauactss and RyRtauinact (not shown), while CaT-tdecay is affected by INaCamax, INaKmax, and GCab (not shown).

The alternans markers also revealed different sensitivities in the nSR and cAF populations. The alternans threshold and range in the nSR population are mostly related to GCaL, INaKmax, GCab, and RyRtauinactss, while in the cAF population they are mostly related to GCab, GK1, GKur, and INaKmax. ΔAPD_max, on the other hand, depends mostly on GCaL GK1, GCab, and GKur (not shown).

Figure 8 shows the global sensitivities of the AP, CaT, and alternans markers as percentages. This result highlights the relative importance of the various ionic currents on the overall AP morphology, with GNa, GK1, GKur, and GCab having the most influence on AP properties. Both CaT morphology and alternans behavior displayed overall greatest sensitivity to GCaL, GCab, INaKmax, and INaCamax. Our results also evince differences in sensitivities of the nSR and cAF populations towards some of the tested parameters, in particular towards INaKmax, INaCamax, GCaL, and GCab.

The nSR_tuned baseline obtained from fine-tuning model parameters displays no APD alternans under dynamic pacing, as designed. AP traces obtained with steady-state pacing at 1000 ms PCL show similar AP morphology as the original nSR baseline, with a slight APD90 prolongation of 35 ms. This prolongation in APD90 is a consequence of the models selected into the new subpopulation all displaying a longer APD90 than the original baseline. However, it should be noted that the scheme chosen in this study for the modulation of parameters to produce a model variant with no alternans is not unique. Other parameter fine-tuning approaches we tested could also produce model variants with no APD alternans under dynamic pacing, these showing no prolongation and even slight shortening of APD90. This attests that model parameters can be fine-tuned within certain tight ranges to produce desired outcomes. However, this observation also seems to suggest that there exists no direct link between APD90 and suppression of APD alternans during dynamic pacing. This observation can be related to the fact that the model parameters changed most dramatically to produce model nSR_tuned, RyRtauinact, iNaKmax, and iCaPmax, and are more closely related to the calcium cycling system than to the AP duration. It can thus be inferred that APD alternans underpinned by Ca-cycling abnormalities will not be suppressed by modulation of APD in this context.

The nSR population presented greater variability in AP and CaT morphology as compared to the cAF population. In particular, APD90 showed the greatest difference in variability between nSR and cAF. The cause might be associated with the electrical remodeling in cAF, such as the 68% increase in I_K1. In a previous study,²⁹ it was found that the hyperpolarization induced by the increased I_K1 in cAF reduces cell excitability, resulting in the stabilization of RMP. Moreover, the authors found that the cAF model was less sensitive to changes in the model parameters than the nSR model, a result which is confirmed by our results.

Our SA showed that AP morphology was mostly affected by GK1, GNa, GCab, Gto, GKur, GCaL, and INaKmax, which are all conductances which scale major cell membrane currents. The collective importance of these major depolarizing and repolarizing currents in the modulation of AP morphology, in particular the five currents I_K1, I_to, I_CaL, I_Kur, and I_NaK, has also been shown by other authors.^16,22,30,31

Increased calcium (depolarizing) and decreased potassium (repolarizing) currents caused APD prolongation in both nSR and cAF, Fig. 7, indicating the dominant role of depolarization in extending the AP. Increased GNa, however, resulted in APD shortening, which may be attributed to early activation of repolarizing currents. Increased GNa also resulted in increased APA and dVdt_max, as expected (I_Na is the main depolarizing current of the AP). RMP showed a negative correlation with GK1; I_K1 is a repolarizing current and, therefore, when it increases, the membrane becomes more hyperpolarized. A more negative RMP allows more sodium channels to recover from inactivation earlier, and to be available at the onset of the next beat,²³ resulting in a higher upstroke velocity. This is reflected by the positive correlation between dVdt_max and I_K1.

In our analyses, larger I_CaL increased ΔCaT and prolonged APD. I_CaL stimulates Ca²⁺ release from the sarcoplasmic reticulum (SR) and is deactivated by the rise in intracellular Ca²⁺. Therefore, a much larger ΔCaT could result in early inactivation of I_CaL. Additionally, increased Ca²⁺ content in the cytosol can activate I_NCX, which pumps Na⁺ inside the cell and activates I_NaK, promoting repolarization and APD prolongation.

Our analyses also showed that an increase in I_NaK is correlated with larger APA and dVdt_max, possibly through intracellular Na⁺ homeostasis. I_NaK has earlier been shown to be linked to calcium cycling and rate adaptation in a model of the atrial cell.¹⁵ I_NaK and I_NCX have additionally been reported to be associated with altered calcium cycling.^15,29 In these analyses, CaT₀ and CaT-tdecay decreased with increasing INaCamax, which reflects the normal function of NCX of extruding calcium from the cell. Decreases in both INaCamax and INaK also caused an increase in Alt-th (i.e., alternans occur at slower pacing rates) and Alt-range, with INaCamax having a more pronounced effect and INaKmax a less pronounced effect on nSR as compared to the cAF population (Fig. 8). This highlights the importance of I_NCX and I_NaK in rate adaptation, in nSR and cAF populations, as previously reported.²⁹ 

In the nSR population, an increase in the inactivation time constant of the RyRs, implying less time occupying an open state, also resulted in increased Alt-th (increase in the PCL at which alternans occurred), Alt-range (range of PCL at which alternans occur), and ΔAPD_max (difference between short and long beats when alternans are present as a measure of degree of alternans). RyRs, as the primary intracellular Ca release channel, is one of the main components implicated in calcium cycling and alternans behavior: Chang et al.³² found that decreasing the inactivation rate constant of RyR resulted in the promotion of APD alternans. In our analyses, we found that an increase in the inactivated state of the RyRs may promote alternans, possibly by allowing more calcium to be released from the SR into the cytosol and prolonging the CaT. If calcium accumulates in the cytosol, the SR may not be fully repleted before the next beat, resulting in decreased CaT. A smaller amount of calcium in the cytosol, in turn, allows enough time for calcium to be re-uptaken into the SR by the SERCA pumps, thus restarting the cycle and leading to CaT alternans. As a result, APD alternans can occur through calcium-voltage coupling. This effect was not observed in the cAF population possibly due to the reduced I_CaL preventing cytosolic calcium overload. However, a more in-depth analysis is required to further investigate these observations.

GCab also seems to have had a greater effect on AP properties than expected, as this is a background current not directly based on experimental measurements. The high sensitivity of the model to GCab across all biomarkers is somewhat surprising and requires further investigation. It is possible that the large dependency on GCab is an artifact resultant of an overestimation of the magnitude of I_Cab in the model.

Note the increased sensitivity of the cAF population to GNa as compared to the nSR population. This is likely related to increased GK1 in cAF which results in a slight hyperpolarization of the RMP in the cAF population. This slightly more negative RMP is sufficient to have a great impact on sodium channel availability, as demonstrated by Skibsbye et al.²³ 

The two simulated populations also displayed distinct dynamic behavior. The cAF population revealed enhanced stability of the restitution curves, showing the steady tendency for the APD90 to decrease as pacing frequency increased (PCL decreased), whereas the restitution curves of the nSR population showed more erratic behavior with increased pacing frequency. This was somewhat unexpected given the higher instability and impaired rate adaptation in cAF in clinical reports.^25–27 Note, however, that all referenced studies were conducted in whole hearts and not on single cells. It has been shown that the electrical restitution properties can change substantially in the presence of electrotonic effects (result of cell-cell coupling) and cardiac memory. Thus, it is expected that alternans behavior at the single cell level differs from that of the whole tissue, where cell-cell coupling can alter or even suppress alternans.^33,34 Alternans at the tissue level are induced by additional effects such as structural remodeling,³⁵ dispersion of refractoriness, and conduction velocity restitution,^7,34 which are not directly related to cell membrane stability. Therefore, these results lend weight to the hypothesis that the electrical remodeling associated with cAF at the cell level, which accounts only for alterations in membrane channels and calcium cycling transporters, may be “cardioprotective” in the sense that these alterations stabilize the cell membrane and thereby reduce alternans and dynamical instability at the level of the single cell.

On the other hand, due to cardiac memory, steady state is highly dependent on the pacing duration and on the initial state of the system. In order to make a direct comparison between different models, the pacing duration at each PCL during the dynamic restitution protocol was kept constant. The consequence of using this protocol is that it is not possible to distinguish between steady state alternans and transient alternans. That is, this dynamic restitution experiment detects alternans, regardless of whether these alternans are transient or stable. Thus, a possible extension of the method could involve establishing a dynamic restitution protocol with an adaptive pacing duration to guarantee that only persistent alternans are included in the analysis.

In the present study, we refined a human atrial single cell model by selecting a subpopulation of models with an alternans threshold at higher pacing rates and longer APD90 from an nSR population, which resulted in a new, more physiologically relevant ‘average’ baseline model with no APD alternans. We also generated two populations of models based on the Koivumäki human atrial cell model,²³ nSR and cAF, by varying parameters related to AP morphology and calcium cycling, such as maximum ion channel conductances, pump fluxes, RyR gating variables, and SR calcium leak. A series of in silico experiments were then performed in order to measure AP and CaT markers, as well as dynamic restitution curves of each model. The populations were functionally calibrated to exclude outliers and several cell biomarkers related to AP and CaT morphology and APD alternans behavior were measured in order to characterize cell stability. The functionally calibrated populations were also subject to SA in order to quantify the sensitivity of different biomarkers to the varied parameters.

Sensitivity analysis with respect to 22 parameters showed that the nSR and cAF populations have high sensitivities with respect to similar parameters related to the calcium currents, I_CaL and I_Cab, I_K1, I_Na, I_Kur, I_NaK, and I_NCX. We have also found that the nSR population exhibited a larger variability of APD90 and dynamic restitution behavior than the cAF population, likely due to the stabilizing effect of electrical remodeling in cAF, as outlined in Secs. IV B and IV D.

The method presented here could be extended to include other single-cell and tissue arrhythmia biomarkers, such as EADs/DADs, effective refractory period, conduction velocity restitution, tissue-level alternans, and spiral wave inducibility. These phenomena have also been linked to arrhythmogenicity and atrial fibrillation initiation and maintenance, and therefore constitute logical extensions of the method. Furthermore, future work will include the definition of a dynamic stability score based on a comprehensive set of arrhythmia biomarkers to rate or classify models according to their stability and propensity to display arrhythmogenic behaviors. We have also observed that the alternans threshold exhibits an interdependence with other cell biomarkers, resulting in two distinct clusters for the nSR and cAF populations, as shown in Fig. 9. This suggests that Alt-th is a rather complex mechanism that results from an interplay of multiple cell behaviors, which in turn exhibit some level of interdependence among them. This highlights that more advanced classification via non-linear regression methods may provide additional insight into arrhythmogenic cellular mechanisms. The prospective dynamic stability score could also be combined with cluster analysis to classify or segment between populations or subpopulations with distinct biomarker targets. The resulting clustered population biomarkers would then represent, for instance, “arrhythmogenic” or “non-arrhythmogenic” models, thus generating a threshold score between the two clusters, which would allow us to better classify models according to, e.g., arrhythmia risk based on specific known cell prescribed biomarkers.

This project has received funding from the European Union's Horizon 2020 MSCA ITN under Grant Agreement No. 675351. We wish to thank the founder of the Allen Institute for Cell Science, Paul G. Allen, for his vision, encouragement and support. We also wish to thank Valeriya Naumova for critical reading of the manuscript.