A population of coupled whole-ventricle and torso models was built taking as a basis the computed tomography (CT) data of a patient (Kania et al., 2017), as shown in Figure 1. The model included the ventricular wall, ventricular and atrial cavities, torso surface, lungs and an approximate anisotropic skeletal muscle layer, which were segmented from the CT data. A hexahedral mesh of the heart with 200 μm resolution and a torso mesh with 1 mm resolution were created.

Ventricular electrical propagation was simulated with a monodomain reaction-diffusion model: ∂ V m ∂ t = 1 C m β − 1 ∇ ⋅ G ∇ V m − I ion , (1) where V _m is the transmembrane potential, C _m the membrane capacitance, β the membrane surface-to-volume ratio (the amount of membrane found in a given volume of tissue), I _ion the sum of all transmembrane ionic currents, and G the monodomain conductivity tensor G _int G _ext/(G _int + G _ext) with G _int and G _ext representing the intracellular and extracellular conductivity tensor fields, respectively. The fiber orientations used to compute the conductivity tensors were assigned with a rule-based method (Streeter et al., 1969; Potse et al., 2006). β was set to 800 cm⁻¹ for myocardium (Potse, 2018). The G _int and G _ext tensors were defined according to the following values: intracellular transverse and longitudinal conductivities were 0.3 and 3 mS/cm and extracellular transverse and longitudinal conductivities were 1.2 and 3 mS/cm, respectively (Potse, 2018).

Temporal integration was done with a forward Euler scheme with a time step of 0.01 ms and a spatial resolution of 0.02 cm. For accuracy, gating variables in the membrane model were integrated with the Rush-Larsen method (Rush and Larsen, 1978), with a time step of 0.01 ms. The simulations were performed using a recent version of the Propag-5 software (Krause et al., 2012) at a stimulation frequency of 1 Hz.

Cellular electrophysiology was represented by the human ventricular myocyte model of Ten Tusscher and Panfilov (2006). The updates to the Ten Tusscher-Panfilov model published by Severi et al. (2009a) were incorporated to adequately represent the relationship between AP duration (APD) and [Ca²⁺]. The L-type calcium current was described as: I CaL = G CaL d f f 2 f Cass 4 V m − 15 F 2 R T × [ C a 2 + ] − 0.25 [ C a 2 + ] ss e 2 V m − 15 F / R T 1 − e 2 V m − 15 F / R T , (2) where G _CaL is the maximal I _CaL conductance, d is a voltage-dependent activation gate, f and f ₂ are voltage-dependent inactivation gates, f _Cass is an inactivation gate dependent on the calcium concentration in the dyadic space (in mM), denoted by [ C a 2 + ] ss , F is the Faraday constant, R is the gas constant, T is the temperature and V _m is the transmembrane potential (in mV). The updated formulations for f ₂ and f _Cass of Severi et al. (2009a) read d f Cass d t = k f Cass,inf − f Cass τ f Cass , (3) where k = 0 if f _Cass,inf > f _Cass and V > − 60 mV, and k = 1 otherwise, f Cass,inf = 0.9 1 + exp [ C a 2 + ] ss − 1.95 0.15 + 0.1 , (4) τ f Cass = 80 1 + [ C a 2 + ] ss 0.05 2 + 1 , (5) and f 2 , i n f = 0.3 1 + e V m + 35 7 + 0.7 , (6) τ f 2 = 600 exp − V m + 25 2 170 + 31 1 + exp 25 − V m 10 + 16 1 + exp V m + 30 10 . (7)

The initial state for each simulation was pre-calculated from a single cell simulation, one for each cell type: endocardial, midmyocardial (M cells) and epicardial. The values of the model state variables after 1000 paced beats were considered as representative of the cell at steady state.

A total of 7 whole-ventricle models with variations in the proportions of endocardial, midmyocardial and epicardial cells were simulated, with the thickness of endocardial and midmyocardial layers ranging from 10% to 50% and the epicardial layer from 20% to 60%. The different cell types differ in the ion channel characteristics, as described in Ten Tusscher et al. (2004) and Ten Tusscher and Panfilov (2006). We used the notation Cuvw, where C stands for the word “case” and u, v and w denote the thicknesses of the endocardial, midmyocardial, and epicardial layers in tenths of the total wall thickness, respectively (e.g., C532 represents the case with 50%, 30%, and 20% of endocardial, midmyocardial and epicardial cells, respectively) (Bukhari et al., 2021). The population of ventricular models included the following combinations of transmural heterogeneities: C136, C154, C316, C334, C352, C514, C532.

The extracellular potential, ϕ _ext, was computed by solving ∇ ⋅ G int + G ext ∇ ϕ ext = − ∇ ⋅ G int ∇ V m (8) in the torso model, with V _m simulated by the monodomain reaction-diffusion model (1) (Potse, 2018). Since we needed to know ϕ _ext only at a few locations for the computation of the ECG we used a Green’s function of the operator ∇ ⋅ ((G _int + G _ext)∇.) for each of these locations to solve this equation efficiently. The ECG leads can be represented by a linear combination of Green’s functions because each recorded ECG lead is a linear combination of ϕ _ext at two or more points. These linear combinations of Green’s functions are termed lead fields (McFee and Johnston, 1954; Colli-Franzone et al., 2000). Our lead fields were computed by solving an equation similar to (8) but as a lead field needs to be computed only once for each ECG lead this approach is much more efficient than solving (8) for each time step in each simulation (Potse, 2018). Five-beat ECGs were simulated at a sampling frequency of 1 kHz.

Principal component (PC) analysis was spatially applied to the T waves (or QRS complexes) of the eight independent leads (Castells et al., 2007) to enhance the T wave (or QRS complex) energy. For each model in the population, the coefficients defining the PC transformation were obtained from the eigenvectors of the 8 × 8 inter-lead auto-correlation matrix estimated by including all segmented T waves (or QRS complexes) at the physiological values of [K⁺], i.e., 5.4 mM, and [Ca²⁺], i.e., 2 mM. Individual ECG recordings of the same model corresponding to other values of [K⁺] and/or [Ca²⁺] were then projected onto the direction of the first PC determined as described above and used for further analysis. We chose to analyze the transformed lead, which represented the main characteristics from the combination of all leads and allowed for an improved morphological characterization.

A wavelet-based single-lead delineation method was used for QRS detection and wave delineation of each of the 12 leads (Martínez et al., 2004). Similarly, the T waves and QRS complexes in the first PC were delineated using the single-lead delineation algorithm described by Martínez et al. (2004). The onset, peak and end of the T waves and QRS complexes were determined and used for subsequent computation of T wave and QRS complex markers.

Simulation results were compared with the results obtained by applying the same analysis to ECG recordings from patients. The study population included 29 ESRD patients from Hospital Clínico Universitario de Zaragoza, Spain, from which 48-h 12-lead ECGs were acquired (see Supplementary Table S1 in the Supplementary Material with [K⁺], [Ca²⁺] and RR values in the patients). T waves (QRS complexes, respectively) in the first PC were obtained and delineated to compute the below described markers (Palmieri et al., 2021a; Bukhari et al., 2022a; Bukhari et al., 2022b).

To assess the relationship between [K⁺], [Ca²⁺] and their combination with T wave and QRS complex characteristics in each model of the population, simulations were conducted under varying values of the electrolytes. The range of simulated [K⁺] values included the default level in the Ten Tusscher-Panfilov model, i.e., [K⁺] = 5.4 mM, as well as other levels below and above it: [K⁺] ∈ {3, 4, 5.4, 6.2} mM. In the case of [Ca²⁺], the range of simulated values included the default level of 2 mM and values around it: [Ca²⁺] ∈ {1.4, 2, 2.6, 3.2} mM. The combinations of [K⁺] and [Ca²⁺] included: (3, 3.2), (4, 2.6), (5.4, 2.0), and (6.2, 1.4) mM. The simulated ranges are similar to those observed in patients during HD (Palmieri et al., 2021a; Bukhari et al., 2021; Bukhari et al., 2022b).

Linear Pearson correlation analysis was also performed to assess the effects of [K⁺] and [Ca²⁺] on each investigated T wave and QRS complex marker.

Sensitivity analysis was performed to quantify how the proportion of endocardial, midmyocardial and epicardial cell layers modulated inter-individual variability in simulated T wave and QRS complex morphology markers at different [K⁺], [Ca²⁺] or their combination levels.

For each T wave and QRS complex marker at each given concentration of [K⁺] ([Ca²⁺] or both [K⁺] and [Ca²⁺], respectively), the percentage of change D Y ; c ; a i in marker Y and its sensitivity S _Y;c;a1,a2 to changes in the proportion of cells was computed as follows (Romero et al., 2009; Bukhari et al., 2019; Bukhari et al., 2021): D Y ; c ; a i = Y c ; a i − Y control Y control ⋅ 100 , i ∈ 1,2 (9) S Y ; c ; a 1 , a 2 = D Y ; c ; a 2 − D Y ; c ; a 1 100 a 2 − a 1 , (10) where Y _c;a is computed by averaging the changes for all possibles combinations having a percentage a _i of cell layer c, c ∈ {endo, mid, epi}, with respect to the reference taken from case C532 (Y _control). The values of a ₁ and a ₂ were taken as the minimum and maximum proportions of cells in each layer, respectively: 10% and 50% for endocardial and midmyocardial cells, and 20% and 60% for epicardial cells.

T wave markers (T _w, T _S/A, d w , T u , d _a,T, d w , T NL , d a , T NL ) computed from simulated ECGs at varying [K⁺], [Ca²⁺] and their combinations are shown in Figure 4. The morphological markers d w , T u , d _a,T, d w , T NL and d a , T NL changed markedly at varying [K⁺] and [Ca²⁺]. Large differences between models in the population could be observed for all analyzed T wave markers. Supplementary Figure S2 in the Supplementary Material shows changes in QRS complex markers at varying [K⁺] and [Ca²⁺] and their combinations, with high variability between models for all markers.

Figure 5 shows T waves and the analyzed markers T _w, T _S/A, d w , T u , d _a,T, d w , T NL and d a , T NL , computed from the ECGs of a particular model, C514, and a particular patient, P10, when concomitantly varying [K⁺] and [Ca²⁺]. More peaked T waves can be observed with increasing [K⁺] and decreasing [Ca²⁺] in both the model and the patient. Analogously, Supplementary Figure S3 in the Supplementary Material shows QRS complexes and analyzed QRS markers for a simulated case and a patient for simultaneous variations in [K⁺] and [Ca²⁺].

Figure 6 shows a comparison of the changes in the marker d w , T NL when varying both [K⁺] and [Ca²⁺] in the simulated cases and the patients. As can be observed by comparing panels a and b, panels d and e and panels g and h, the models in the population reproduced some specific patterns of change of d w , T NL in the patients, albeit with some quantitative differences. These results were confirmed by computation of correlation coefficients, as shown in panels c, f and i of the same figure. The ability of our in silico population to reproduce d w , T NL trends measured in some of the patients was equally valid for other T wave markers even if they did not present as remarkable changes as d w , T NL when varying [K⁺] and [Ca²⁺]. Analogous results in simulated and patients’ ECGs are depicted in Supplementary Figure S4 of the Supplementary Material for the QRS-based d a , Q NL , which was the one showing the largest changes in response to electrolyte variations.

To assess the extent to which our population of models could reproduce the inter-patient variability in T wave markers at concomitantly varying electrolyte levels, a correlation analysis was performed. Figure 7 shows the Pearson correlation coefficient r between each T wave marker and [K⁺], or [Ca²⁺], in the simulated and the patients’ ECGs. T _w, d w , T u and d a , T NL were the markers most strongly correlated with [K⁺] (median r being −0.70, 0.87, 0.91 in simulations and −0.92, 0.93, 0.75 in patients, respectively) and [Ca²⁺] (median r being 0.70, −0.85, −0.91 in simulations and 0.79, −0.84, −0.75 in patients, respectively). Inter-individual variability in the correlation coefficients associated with T _w and T _S/A was high in both simulations and patients. For all other T wave morphology markers, the variability between models only partly reproduced the variability between patients. Table 1 provides the results for the quantitative comparison between simulated and patients’ T wave markers, in terms of median and interquartile range of r with [K⁺] and [Ca²⁺]. As can be seen from the table, all the analyzed morphology-based T wave markers correlated strongly with [K⁺] and [Ca²⁺] in simulations and patients, with part of the inter-patient variability being reproduced by the models.

Supplementary Figure S5 and Supplementary Table S2 in the Supplementary Material show correlation coefficients between QRS complex markers and electrolyte levels in simulated and patients’ ECGs.

The results of the sensitivity analysis performed to investigate how different proportions of endocardial, midmyocardial and epicardial cells may explain individual T wave responses when varying both [K⁺] and [Ca²⁺] are presented in Table 2 for all the analyzed simulated T wave markers. The highest sensitivity values were observed for morphology-based T wave markers, particularly d w , T u and d w , T NL , when varying the proportion of epicardial cells. Supplementary Table S2 in the Supplementary Material shows sensitivity results for QRS complex markers.

We measured commonly used QRS complex and T wave markers describing characteristics related to their duration and amplitude, including T wave and QRS complex widths (T _w, QRS _w), QRS complex amplitude (QRS _a) and T wave slope-to-amplitude ratio (T _S/A), as well as our recently proposed morphological variability-based markers ( d w , T u , d w , Q u , d _a,T, d _a,Q, d w , T NL , d w , Q NL , d a , T NL , d a , Q NL ) (Ramírez et al., 2017; Palmieri et al., 2021a; Bukhari et al., 2021; Bukhari et al., 2022a). We evaluated these markers in simulated ECGs derived from coupled heart-torso models with different proportions of endocardial, midmyocardial and epicardial cells, at varying [K⁺] (3–6.2 mM), [Ca²⁺] (1.4–3.2 mM) and their combinations.

We simulated a range of [Ca²⁺] values to assess the impact on the ECG markers. In the patients, total serum calcium concentration ([Ca²⁺]) was measured, which has higher values than the corresponding ionized calcium concentration. The values of [Ca²⁺] in our analyzed ESRD patients were in the ranges reported in other studies (Siyam and Klachko, 2013; Janmaat et al., 2018; Noordam et al., 2019; Akter et al., 2022). Supplementary Table S1 in the Supplementary Material shows the [K⁺], [Ca²⁺] and RR interval values in the 29 patients. Changes in the RR interval (inverse of heart rate) from the start to the end of the HD session were not statistically significant [see Table 2 in Bukhari et al. (2021) for the effects of RR intervals on all the analyzed T wave markers and Table 2 in Bukhari et al. (2022a) for the effects of RR intervals on all the analyzed QRS markers in the patients]. Consequently, we decided not to explore the effects of other stimulation frequencies in our heart-torso simulations.

We found that most of our simulated QRS and T wave markers (all except for QRS _w, d _a,T and d _a,Q) presented a diversity of patterns in their relationships with [K⁺] and [Ca²⁺] that were in line with our observations in ESRD patients during and after HD (Palmieri et al., 2021a; Bukhari et al., 2021; Bukhari et al., 2022b; Bukhari et al., 2022b). Overall, the T wave morphology markers d w , T u , d w , T NL and d a , T NL were the ones that most notably changed with [K⁺], somewhat less notably with [Ca²⁺] and very remarkably with their combination, in average over the population of in silico models (Figure 4). When we tested lower [Ca²⁺] values of 0.8, 1, and 1.2 mM, which would correspond to ranges of ionized calcium concentrations, the changes in the T wave markers were accentuated with respect to the ones observed for [Ca²⁺] varying from 1.4 to 3.2 mM, as illustrated in Supplementary Figure S6 of the Supplementary Material. Nevertheless, the variations in the T waves when only varying [Ca²⁺] were still smaller than the ones measured when varying [K⁺] and [Ca²⁺] simultaneously. Our simulation results were in line with averaged changes presented by the analyzed markers in ESRD patients, as shown here and in previous studies where the markers were first evaluated during HD (Palmieri et al., 2021a; Bukhari et al., 2021; Bukhari et al., 2022b). This generally corresponded to more peaked T waves, with lower amplitude and/or longer duration, which presented increasingly larger morphological differences with respect to a reference T wave (calculated at physiological electrolyte levels) when [K⁺] was increased and [Ca²⁺] was decreased (Figure 5) (Palmieri et al., 2021a; Bukhari et al., 2021; Bukhari et al., 2022b). Regarding the QRS complex, markers describing amplitude characteristics, like QRS _a and d a , Q NL , also varied appreciably with [K⁺], [Ca²⁺] and, particularly, with their combination, in both the simulations and the patients. Typically, the averaged QRS complex in the simulated and patient populations became lower in amplitude and had larger non-linear amplitude morphological variations than a reference physiological QRS complex when evaluated in response to increasing [K⁺] and decreasing [Ca²⁺] (Supplementary Figure S3 in the Supplementary Material) (Bukhari et al., 2022a).

For confirmation of the ability of our in silico population to describe averaged trends of T wave and QRS complex changes with electrolyte variations in patients, a correlation analysis was performed. Among the analyzed T wave markers, d w , T u and d a , T NL were the most strongly correlated with [K⁺] and [Ca²⁺] in both the simulations and the patients (absolute median correlation coefficient ranging from 0.75 to 0.93). In the case of QRS complex markers, QRS _a, d w , Q NL and d a , Q NL were the ones most strongly associated with [K⁺] and [Ca²⁺] in the simulations and the patients (absolute median correlation coefficient ranging from 0.71 to 0.98). These results support the use of our in silico population of models to describe the average response of T wave and QRS complex to variations in electrolyte levels like those seen in ESRD patients during and after HD.

After assessing how our population of in silico models could reproduce the general trends of electrolyte-induced changes in ECG markers observed in ESRD patients, we investigated how it may explain the high inter-patient variability in such changes. In agreement with the diversity of patterns of QRS complex and T wave changes in response to [K⁺] and [Ca²⁺] variations measured in the patients, we found different ventricular models presenting largely distinct behaviors that covered some of the patterns described in the patients. Models with thin epicardial layers were mainly related with a non-monotonic behavior of the evaluated ECG markers similar to that presented by some ESRD patients. This was particularly clear for the T wave morphology marker d w , T NL , illustrated in Figure 6. Most models with thick epicardial layers reproduced the increasing linear trend of d w , T NL observed in many patients (Figure 6). For the QRS complex, we found that models with thick midmyocardial layers could replicate the increasing linear trend in the marker d a , Q NL measured in some patients (Supplementary Figure S4 of Supplementary Material). Also, the non-monotonic behavior of d a , Q NL in many other patients was obtained for models with thick epicardial layers (Supplementary Figure S4 of the Supplementary Material). Therefore, we may conclude that the different transmural heterogeneities simulated in this study may explain the inter-patient variability in the ECG response to [K⁺] and [Ca²⁺] variations.

To quantify how different ventricular wall compositions may explain individual ECG responses at a range of [K⁺] and [Ca²⁺], we computed QRS and T wave marker sensitivities to variations in the thickness of the endocardial, midmyocardial and epicardial cell layers. The highest sensitivity values were observed for T wave morphology markers d w , T u , d w , T NL and d a , T NL and for QRS morphology markers d w , Q NL , d _a,Q and d a , Q NL . Extremely high sensitivities were found for some of these T wave markers at abnormally high levels of [K⁺] and low levels of [Ca²⁺], which may explain the remarkably large inter-patient variability observed at the beginning of the HD, when [K⁺] was enhanced and [Ca²⁺] was reduced. In particular, d w , T u very notably increased for progressively larger proportions of epicardial cells and it decreased for progressively larger proportions of endocardial or midmyocardial cells. The largest sensitivity of d w , T u to variations in the proportion of epicardial cells within the ventricular wall agrees with previous reports showing the contribution of epicardial cells to other forms of repolarization variability like T wave alternans (Janusek et al., 2018). For [K⁺] = 6.2 mM and [Ca²⁺] = 1.4 mM, sensitivity values above 110% were found, which corresponded to a coefficient of determination of 0.93 for the relationship between d w , T u and the proportion of epicardial cells. Previous experimental and theoretical studies have described how cell type distributions influence ECG characteristics (Hooft van Huysduynen et al., 2005; Okada et al., 2011; Rivolta et al., 2015; Ramírez et al., 2017; Janusek et al., 2018; Srinivasan et al., 2019). In particular, Janusek et al. (2018) reported that the contribution of epicardial cells to T wave alternans was significantly higher than that of midmyocardial cells, which would be in line with our results for T wave morphological variability. In Rivolta et al. (2015), cell type distributions were shown to highly affect both repolarization and T wave morphology parameters, in concordance with our findings on the impact of ventricular wall composition on T wave characteristics, particularly at abnormal [K⁺] and [Ca²⁺] values.

For the above described analysis, we included different proportions of endocardial, midmyocardial and epicardial layers that comprised those used in previous simulation studies, which were in turn based on experimentally reported transmural cell proportions (Drouin et al., 1995; Gima and Rudy, 2002; Ten Tusscher and Panfilov, 2006; Bueno-Orovio et al., 2008; Pueyo et al., 2010; Pueyo et al., 2011). In particular, the following proportions were used for endocardial, midmyocardial and epicardial layers in previous simulation studies: Gima and Rudy (2002) used 4:2:4; Bueno-Orovio et al. (2008) used 5:3:2 and Pueyo et al. (2010) used 3:3:4. In a previous study (Bukhari et al., 2021), we included larger ranges of cell type proportions in simulated 1D fibers, with the endocardial percentage ranging from 10% to 50%, the midmyocardial percentage from 10% to 50% and the epicardial percentage from 20% to 80%. Here, we included a subset of those cases, which covered the proportions most commonly employed in previous simulation studies. Even if we showed that variability in the thickness of endocardial, midmyocardial and epicardial layers may be a cause for the observed inter-individual variability in the ECG response to electrolyte variations, this was based on the large range of cell type proportions simulated in our study. If the variation in the myocardial cell distribution in the patients were substantially smaller, the differences in ventricular wall distribution might play a less relevant role in the observed inter-patient variability of the ECG response to electrolyte variations. Moreover, it is not even known whether M cells with distinct electrophysiological properties exist or not. There has been and there is still controversy regarding these cells, with some studies showing evidence of their existence (Liu et al., 1993; Drouin et al., 1995; Anyukhovsky et al., 1996; Baláti et al., 1998; Stankovicova et al., 2000; Yan et al., 2003; Wilson et al., 2011) and others questioning it, particularly in intact hearts (Rodríguez-Sinovas et al., 1997; Bryant et al., 1998; Conrath et al., 2004; Janse et al., 2012).

Our in silico analysis could pave the path to understand how patients with different ventricular wall compositions may present largely different responses to serum electrolyte variations, even if other factors, such as interventricular, apicobasal or anteroposterior differences (Janse et al., 2012; Opthof et al., 2017), might also explain inter-patient variability. These results could help to improve monitoring serum ion concentrations and predicting arrhythmic events in ESRD patients based on investigations combining in silico modeling and simulation with signal processing of the ECG.

Several previous studies have characterized ECG features in relation to electrolyte concentrations. Hernández Mesa et al. (2018) and Pilia et al. (2017) computed ECGs at different [Ca²⁺] from the ventricular electrophysiology and a torso model. No changes with [Ca²⁺] were found in QRS duration, while the R wave amplitude and energy diminished with decreasing [Ca²⁺] (Hernández Mesa et al., 2018), which would agree with our observations of reduced QRS _a and d a , Q NL . Also, the T wave slopes were reported to increase with decreasing [Ca²⁺] (Hernández Mesa et al., 2018) and increasing [K⁺] (Pilia et al., 2017), which is in line with our observations of more peaked T waves at low [Ca²⁺] and high [K⁺].

Other studies have simulated human ventricular electrophysiology in a one-dimensional transmural model and have derived pseudo-ECGs at varying [K⁺] (Kharche et al., 2012; Corsi et al., 2017). Increases in T wave slope were reported in association with hyperkalemia, which is in agreement with our results.

Inter-patient variability in the relation between electrolytes and ECG markers has been reported by several authors (Buemi et al., 2005; Montague et al., 2008; Severi et al., 2009b; Alabd et al., 2011; Green et al., 2013; Corsi et al., 2017; Yoon et al., 2018; Pilia et al., 2020; Palmieri et al., 2021b; Palmieri et al., 2022). We are not aware of previous work explaining this variability in terms of variations in ventricular wall composition, except for our previous study in simulated ventricular 1D fibers (Bukhari et al., 2021). By extending our investigations from 1D fibers to 3D ventricular and torso models, we could compute ECG features in a more realistic way that allowed for a more accurate comparison with patient data.

Our simulations were performed with the most realistic available techniques to model the electrophysiology of the heart and the resulting ECGs. The electrophysiology was represented by a widely-used and detailed model of the human ventricular myocyte (Ten Tusscher and Panfilov, 2006) that is thoroughly rooted in experimental data. We included the updates to the Ten Tusscher-Panfilov model proposed in Severi et al. (2009a) so that the model reproduced a physiological APD-[Ca²⁺] relationship as described in several in vitro and in vivo studies (Temte and Davis, 1967; Leitch and Brown, 1996; Bai et al., 2005). This model allowed running cell and tissue simulations at an affordable computational time. Initially, we considered the O’Hara-Rudy model (O’Hara et al., 2011) and the Grandi-Pasqualini-Bers model (Grandi et al., 2010), but these models had limitations with respect to the simulation of hyperkalemic effects in tissues, as they did not show propagation of excitation for [K⁺] values larger than 6 and 8 mM, respectively (Dutta et al., 2017). Subsequently, we considered the modified version of the O’Hara-Rudy model proposed by Dutta et al. (2017), but we observed that this model was not able to reproduce the inverse relationship between APD and [Ca²⁺] reported experimentally. Another model recently proposed by Bartolucci et al. (2020) which is a modified version of the O’Hara-Rudy 2011 model, can reproduce a physiological APD-[Ca²⁺] relationship but requires further testing in 3D tissues. The Bartolucci et al. model could be a good option for future studies, considering its ability to reproduce a physiological APD-[Ca²⁺] relationship.

We found that variations in wall composition may explain inter-patient variability in the ECG response to electrolyte changes. While we have shown that variability in the thickness of endocardial, midmyocardial and epicardial layers may be a cause, we cannot exclude additional contributing factors. On top of including other types of ventricular heterogeneities, as discussed above, future studies could also investigate electrophysiological differences across the ventricular wall, e.g., by incorporating transmural differences in the I _CaL formulation. In any case, the findings of our study remain to be confirmed. If confirmed, this could imply that accurate serum ion concentration estimation in individual patients would require a priori knowledge of the individual patient’s cell type distribution. Methods could be developed to infer such distribution from the ECG and blood data acquired from the patient during one or several prior HD sessions. ECG-based estimators of serum ion concentrations could be subsequently derived as a function of the patient’s cell type distribution.

Our in silico population included 7 whole-ventricle models with different cell type distributions across the ventricular wall, which we used to simulate variations in [K⁺] from 3 to 6.2 mM and in [Ca²⁺] from 1.4 to 3.2 mM. Future studies could extend the methods proposed here to add larger numbers of coupled ventricle-torso models and simulate a wider range of [K⁺] and [Ca²⁺].

The effects of other electrolyte concentrations such as magnesium, [Mg²⁺], which has been reported to be possibly involved in observed alterations in the ECG (Noordam et al., 2019; Naksuk et al., 2017; Jhang et al., 2013; van den Bergh et al., 2004), could be accounted for in the simulations if information on the variation of these electrolytes during and after HD were available. In the present study, we did not investigate [Mg²⁺] variations because serum [Mg²⁺] levels were not measured in our patients.

To investigate other sources of inter-individual variability in the ECG response to serum electrolyte variations, future studies could include other ventricular heterogeneities on top of transmural ones, like interventricular, apicobasal or anteroposterior (Janse et al., 2012; Opthof et al., 2017), which might play a relevant role in determining ECG characteristics under normal and abnormal electrolyte concentrations.