Cell to whole organ global sensitivity analysis on a four-chamber heart electromechanics model using Gaussian processes emulators

Cardiac pump function arises from a series of highly orchestrated events across multiple scales. Computational electromechanics can encode these events in physics-constrained models. However, the large number of parameters in these models has made the systematic study of the link between cellular, tissue, and organ scale parameters to whole heart physiology challenging. A patient-specific anatomical heart model, or digital twin, was created. Cellular ionic dynamics and contraction were simulated with the Courtemanche-Land and the ToR-ORd-Land models for the atria and the ventricles, respectively. Whole heart contraction was coupled with the circulatory system, simulated with CircAdapt, while accounting for the effect of the pericardium on cardiac motion. The four-chamber electromechanics framework resulted in 117 parameters of interest. The model was broken into five hierarchical sub-models: tissue electrophysiology, ToR-ORd-Land model, Courtemanche-Land model, passive mechanics and CircAdapt. For each sub-model, we trained Gaussian processes emulators (GPEs) that were then used to perform a global sensitivity analysis (GSA) to retain parameters explaining 90% of the total sensitivity for subsequent analysis. We identified 45 out of 117 parameters that were important for whole heart function. We performed a GSA over these 45 parameters and identified the systemic and pulmonary peripheral resistance as being critical parameters for a wide range of volumetric and hemodynamic cardiac indexes across all four chambers. We have shown that GPEs provide a robust method for mapping between cellular properties and clinical measurements. This could be applied to identify parameters that can be calibrated in patient-specific models or digital twins, and to link cellular function to clinical indexes.


ToR-ORd model
Global sensitivity analysis Table 1 shows the parameters we considered in the analysis on the ToR-ORd model. We included the conductances for all ionic currents: fast and slow sodium currents G N a and G N aL , transient outward potassium G to , rapid and slow delayed potassium rectifier currents G Kr and G Ks , inward potassium rectifier current G K1 , L-type calcium current P Ca , sodium-calcium exchanger G N CX , sodium-potassium pump current G N aK , calcium pump current G Ca , calcium-sensitive chloride current G ClCa , all background currents G N ab , P Cab , G Kb , G Clb ; and the parameters representing calcium handling: calcium release (J rel ) and uptake (J up ) in the sarcoplasmic reticulum, fraction of sodium-calcium exchangers (I N aCa,SS ) and L-type calcium channels (I CaL,SS ) located in the subspace, calmodulin-kinase kinetics parameters (α CaMK , β CaMK , CaM K o ), calcium buffer concentrations ([CMDN], [TRPN], [BSL], [BSR] and [CSQN]) and calcium diffusion rates (τ diff,Ca and τ tr ). For details about the ToR-ORd model, all symbols and default values, the reader is referred to the original publications [1][2][3]. For GPE training and GSA, the range for all parameters was set to ±25% from the default value. The parameter space was sampled with a latin hypercube design with N=2175 points. For each sample, we ran the ToR-ORd model for 100 beats and a cycle length of 1000 ms. The last beat was used to compute the following calcium transient features: 1. diastolic calcium level (Ca diast ) 2. calcium transient amplitude (Ca ampl ) 3. time to peak (TTP) 4. time to reach 90% of calcium transient decay (RT90).
We did not consider transmembrane voltage dynamics in our sensitivity analysis because the activation at the organ level was computed with a reaction-Eikonal model without diffusion. Therefore, the conduction velocities in the Eikonal model are the only determinants of the speed of electrical propagation.
The performance of the trained GPEs is provided in Table 2. For all four calcium features, the average R 2 and ISE scores were >0.98 and >98.0, showing that the GPEs provide accurate predictions for all outputs. A Saltelli sampling with a Sobol base sequence of N base =1000 samples was generated and the GPEs evaluated to compute the total order effects and to identify important parameters for the calcium features. Fig 1A  shows a heatmap of the total effect of the parameters (x-axis) over the output (y-axis). The calcium transient amplitude was mostly affected by the conductance of the L-type calcium channels (P Ca ), which also affects all other calcium features to a lesser extent. The conductance of the sodium-calcium exchanger (G N CX ) had a high effect on the diastolic calcium, which was also affected by the ratio of the L-type channels located in the subspace (I CaL,SS ). The time to peak and the 90% relaxation time were mostly affected by the uptake of calcium in the sarcoplasmic reticulum (J up ). The calmodulin-kinase parameters (α CaMK , β CaMK , CaM K o ) affected all calcium transient features, but with smaller effects. The most important calcium buffer was the troponin ([TRPN]), which impacted the calcium decay. The sarcolemmal binding sites in the subspace [BSL] and the calcium diffusion rate from the subspace to the sarcolemma had a small but significant effect on the time to peak of the calcium, while all other parameters were not significanlty important. Fig 1B shows the ranking of the parameters according to their maximum total effect across all outputs, normalised to sum up to 1. The orange bars represent the most important parameters needed to collectively explain >90% of the output features. All other parameters were excluded from the following analysis, as they did not significantly impact the calcium transient. Table 1. ToR-ORd model parameters. The first three columns show the parameter name, its default value, the GSA and HM ranges and its meaning. The last column provides the original paper the symbol and the default value refer to. The gray rows indicate unimportant parameters as found with the GSA on the ToR-ORd model. Abbreviations and symbols: Na + =sodium, K + =potassium, Ca 2+ =calcium, SR=sarcoplasmic reticulum, CaMK=calcium/calmodulin-dependent protein kinase II, DS=dyadic space.

History matching
In the previous section, we identified the ToR-ORd parameters that mostly affected the simulated calcium transient. In this section, we aim to use HM to find areas in the parameter space where these parameters provide physiological calcium features. The target mean µ and standard deviation σ for the calcium features were based on literature values from datasets from Coppini et al [4] and Piacentino et al [5]. Coppini Table 3, and were computed using the following formulas for mean and standard deviations for combined datasets: where µ 1 , σ 1 , N 1 and µ 2 , σ 2 , N 2 are the mean, standard deviation and number of samples in the two datasets.
In the ToR-ORd model, the diastolic calcium was not calibrated to match any experimental data. Therefore, with default parameters, the simulated diastolic calcium was smaller than the physiological target values (0.074 µM vs 0.1462±0.01661 µM). To achieve physiological diastolic calcium during the HM, we increased the range for P Ca , G NCX and I CaL,SS to ±75% from their default values, as these parameters were the most important for the diastolic calcium ( Fig 1A). The range for the other important parameters was set to ±50%. Table 3 summarises the results for the HM on the ToR-ORd model. The initial GPEs were trained with N=250 Latin hypercube samples, and the first HM wave was run with a 3.5 threshold on the implausibility measure. The subsequent waves were run by enriching the training set for the GPEs with N simul =128 each time. From the fourth wave onwards, the initial dataset was excluded from the training set to improve GPEs performance. From the first to the last wave, the % of plausible points increased from 19.2% to 98.9% (Fig 2B, blue to red areas) and the mean implausibility measure decreased from 4.4 to 2.3, indicating that the predicted values for the output features are on average closer to the target ranges. In addition, GPE uncertainty decreased, as the mean variance ratio between uncertainty prediction and uncertainty on the data got smaller. Fig 2A shows

ToR-ORd and Land model
In the following sections, we aim to identify important parameters for the active tension transient in three different types of simulations: 1) isometric twitch with constant stretch λ = 1.0 (e.g. no strain); 2) isometric twitch with constant stretch λ = 1.1; 3) isotonic twitch based on the Land model equation with cellular passive stress as in Jung et al [6]. We then use HM to find areas in the parameter space where both calcium and isometric tension are physiological and provide physiological dynamics at the whole organ level. All simulations were run for 100 beats with a basic cycle length of 1000 ms.

Global sensitivity analysis
For the GSA on the Land model, we considered the important parameters for the ToR-ORd model ( Table 1, white rows) and all 17 parameters in the Land model. Table 4 summarises the parameters, their meaning, their default values and their range. The default for k u and n T m were taken from [7], where the Land model was coupled to the ToR-ORd model and recalibrated. The default value for all other parameters was based on the original Land model publication [8]. The Land model parameters were constrained between ±50% from their default value, while the ToR-ORd parameters were sampled from the non-implausible region of the last HM wave (Fig 2B,  6. rest active tension (T rest ).
In addition, for isotonic twich simulations, we computed the % of λ shortening (∆λ) to investigate the effect of calcium and tension parameters on the extent of cell contraction.
There is a high variability in isometric active tension values reported in the literature due to different tissue preparations, temperatures and species [9]. Reported values for the peak in active tension range from 10 kPa in mice trabeculae [10] up to 121±35 kPa [11] in the rat trabeculae and 108±13.8 kPa in the cat right ventricular trabeculae. Based on these literature values and to prevent the GPEs to be trained on unphysiological tension transients, we excluded samples that resulted in T peak <10 kPa. In previous models, T rest was constrained to be <1 kPa [12]. We relaxed this constraint to account for the effect of length dependence and excluded samples that resulted in T rest >2 kPa. This provided N=1584, N=1017 and N=1100 samples for isometric twitches with λ = 1.0 and λ = 1.1, and an isotonic twitch, respectively. GPE performance is summarized in Table 5. For all active tension features, the GPEs performed well, with average R 2 and ISE scores above 0.84 and 97.0, indicating that the GPEs provide accurate prediction of model outputs. For the rest tension during an isometric twitch with λ = 1.1, the R 2 was 0.7843, which is lower than those achieved for the other outputs. This was due to highly non-linear behaviour of the Land model for non-zero strains.
The trained GPEs were used to run a GSA on the active tension transient features. To make sure that the GPEs were not evaluated outside their training space, the base sequence for the Saltelli sampling was screened with the GPEs trained during the HM on the ToR-ORd model to exclude samples providing non-physiological calcium. The base sequence sampling was repeated until the base sequence had N sobol >1000 samples. Fig  3 shows the results of the GSA on the ToR-ORd-Land model. For all three simulated scenarios (isometric λ = 1.0 and λ =1.1, and isotonic) the conductance of the L-type calcium channels (P Ca ) and the conductance of the sodium-calcium exchanger (G NCX ) affected the peak in tension, and the minimum and the maximum tension derivatives. The minimum tension derivative was also affected, but to a lesser extent, by the maximum troponin concentration ([TRPN]). The reference tension T ref also had an impact on the peak in tension and its derivatives, while the reference calcium sensitivity ca 50 has a significant effect on all outputs. The rest tension was highly affected by n TRPN and n Tm , as they define the cooperativity between the calcium and troponin, and the calcium troponin complex with unbound cross-bridges. Cross-bridges kinetics parameters k u , µ, T RP N 50 , r s and r w have a significant effect on the active tension time to peak and transient duration, although their relative contribution changes during an isometric vs isotonic twitches. Finally, the velocity dependence parameter A eff affects the active tension rising time and maximum tension derivative only during an isotonic twich, which is expected as during isometric twitches the velocity of contraction is 0. Fig 4A-C shows the parameter ranking according to their maximum total effect across all outputs for the three different simulated scenarios, while Fig 4D ranks the parameters according to their maximum total effect across all three simulation types. The maximum total effects were then normalised to sum up to 1, e.g. 100%. We considered important those parameters that were needed to explain >90% of output variance in all three simulations (e.g. parameters on the left of the dashed line in Fig 4D). All other parameters were excluded from the following analysis.  Table 5. GPEs performance ToR-ORd-Land. R 2 score and ISE for every split of a 5-fold cross-validation, reported for each output for an isometric twitch with λ = 1.0, λ = 1.1 and an isotonic twitch.

History matching
In the section above, we identified calcium and active tension parameters that are important for the active tension transient. In this section, we want to find areas in the parameter space where these provide both physiological calcium and active tension for an isometric twitch with λ = 1.0, and also result in physiological behaviour at the organ level. The top section of Table 6 shows the target mean values µ and standard deviations σ for the active tension features we included in the HM. The target values for the calcium features were set as explained above in the ToR-ORd model HM. The mean target peak in isometric tension was set to 160 kPa, in range with measurements in ventricular myocytes in mammals [11], with a standard deviation of 15 kPa. The duration of the active tension transient was set to 500±20 ms, consistent with pressure transient duration available from the clinical data. Finally, the rest tension target values were set to 0.5±0.1667 kPa which, together with the three-sigma rule, constrains the rest tension to be below 1 kPa, consistently with previous modelling studies [13]. The default values and parameter space bound for all parameters is provided in Table 4 (white rows). The default values for the calcium parameters was set to the center of the non-implausible region from the last HM wave on the ToR-ORd model, while all other parameters were fixed at their default value listed in Table 4. The reference tension T ref was limited between 100 kPa and 200 kPa, to achieve the target value for the peak in tension of 160 kPa. The bottom section of Table 6 summarises the results for the HM iterations. The initial GPEs were trained on 260 simulations, and the first wave was run with a threshold on the non-implausibility measure of 3.5. The threshold was then decreased to 3.0 for all other waves, and the GPE traning set was enriched with N simul =128 additional simulations. To improve GPEs accuracy, we excluded the initial set of simulations from the GPE training dataset for waves 4 and 5. To improve accuracy even further and avoid training the GPEs with unphysiological samples, waves 6 to 9 were run only keeping the training datasets from the last three waves. During this procedure, the percentage on non-implausible points increased from 2% to 99.2%, and the mean implausibility measure decreased from 7.8 to 2.5. Starting from very uncertain GPEs with a mean variance ratio of 111.56, indicating that the GPEs are more than a 100 times more uncertain than the experimental data, the variance ratio decreased to 0.09. By restricting the parameter space to the non-implausible area (Fig 5B, blue to red), the calcium and active tension transients used to train the GPEs go from unphysiological (Fig 5A, blue) to physiological (red).
Since the samples from the non-implausible area of the last HM wave on the ToR-ORd-Land model will be used to run whole organ simulations, we need to ensure that the model behaves physiologically for all samples. To ensure this, we ran one last HM wave discarding the GPEs (e.g. using cell model evaluations rather than GPE predictions and setting the uncertainty of the prediction to 0). We therefore ran 200,000 test samples with the cell model extracted from the non-implausible region of wave 9, evaluated the implausibility measure for all samples, and excluded those which did not provide physiological calcium and tension outputs. Initial tests in the whole heart simulations allowed us to identify two additional constraints that improved numerical stability of the fully coupled simulations. Specifically, we ran the following checks on each sample before running whole organ simulations: 1. isometric twitch λ = 1.1 and basic cycle length of 1000 ms: T rest <2 kPa to ensure that, even for non-zero strains, the rest tension was not too high B A Simulated model dynamics