© The Authors, published by EDP Sciences, 2020

1 Introduction

Ageing is one of the dominant risk factors for cardiovascular diseases. A large body of experimental data has been gathered on cellular remodeling in the ageing myocardium from animals but very few experimental data are available on the age-related changes in the human cardiomyocyte [1, 2]. Mathematical models could be useful for predicting the adverse consequences of age-related changes in the myocardial function.

Recently, populations of human ventricular in silico action potential (AP) models were used to predict the clinical drug-induced arrhythmic risk of Torsade de Pointes (TdP) arrhythmia [3, 4]. The occurrence of cellular repolarization abnormalities (RA) and increasing of AP duration was evaluated in model populations and showed high accuracy (close to 90%) in predicting TdP risk for a number of references compounds [3, 4]. A similar approach can be used for assessing the pro-arrhythmic consequences of various interventions affecting the cellular mechanisms of excitation-contraction coupling. In this paper, we apply the experimentally calibrated population-of-models methodology [3-6] to investigate variability in the response of ventricular cardiomyocytes to age-related changes in the parameters of ionic currents reported in literature [1, 2]. Besides, the population of models based on a human electro-mechanical cellular model developed by our research group [7] is used here for predicting the pro-arrhythmogenic effects of age-related cellular remodelling.

2 Methods

2.1 Control population of models

In this study we used a new combined electro-mechanical model (TNNP+M) of the human ventricular cardiomyocyte [7] based on the TNNP06 ionic model [8] of action potential (AP) and Ekaterinburg model of the mechanical activity and calcium handling in ventricular cardiomyocytes [9]. The ionic model of AP contains a biophysically detailed description of ionic channels, pumps, and exchange currents, detailed equations for the concentration of intracellular sodium, calcium, and potassium. The mechanical part of the model contains equations for myocardial stress, changes in the length of sarcomeres, and the whole cardiomyocyte, including the kinetics of force-generating cross-bridges (Xb) and regulatory proteins troponin C (TnC) providing for mechano-Ca²⁺ and mechano-electric feedback in the cells. For a full description of model, we refer to [7].

Using a methodology developed by B. Rodriguez group [3-6], we constructed a control population of human ventricular electro-mechanical models based on the TNNP+M model with the reference set of parameters. Similarly to [3], we randomly varied 9 model parameters: 8 parameters of the membrane conductance for main transmembrane ionic currents: maximal conductance for the fast Na⁺ current (GNa), transient outward K⁺ current (Gto), rapid and slow delayed rectifier K⁺ currents (GKr and GKs), inward rectifier K⁺ current (GK1), the Ltype Ca²⁺ current (GCaL); and the Na⁺-K⁺ (NKX) and Na⁺-Ca²⁺ (NCX) exchangers (P_NaK, K_NaCa). Additionally, the maximal velocity of the SERCA pump (Vmax_up) was chosen to be varied as it sufficiently affects CT and subsequently mechanical activity in the cells [10]. At the same time, we did not vary any specific mechanical parameters in this study, focusing on possible effects of the electrical parameter variation on the variability in the cellular CT and mechanics.

We randomly generated 20,000 parameter sets using the Latin hypercube sampling scheme with each parameter ranged from 0 to 200% of a reference value used in the baseline model. For each model sample, 200 cycles at a pacing rate of 1 Hz were calculated to achieve a steady state.

To select non-implausible (non-ImP) models for a control population, the calculated models were filtered against experimental biomarkers derived from the AP and CT recordings in isolated cells from undiseased human hearts as referred in [4, 11, 12] (see Table 1). Corresponding characteristics of simulated signals were compared with experimental data, and the models with outputs out of the physiological diapason were excluded from the control population as implausible.

The “control” population of virtual ventricular cardiomyocytes consisted of 240 model samples with biomarkers within the normal physiological range. Signals of AP, CT and isometric force twitch (FT) developed by every model in the population is shown together with distributions for model parameters and biomarkers in Figure 1.

All the simulations presented in this study were performed using Python3.6 language. Differential equations were solved using CVODE solver – a part of the open-source SUNDIALS suite [13] with just-in-time JIT Numba compilation to speedup calculation.

Fig. 1. A control population of human ventricular electro-mechanical models. A Traces of action potential (AP), Ca2+ transient (CT) and force twitch (FT) in the 240 models in the calibrated population. B Ionic parameter profile of the calibrated population, shown as scaling factors with respect to the baseline model. Central mark at box-plots corresponds to the median, box limits are the 25th and 75th percentiles, and whiskers extend to the most extreme data points. C Distributions of AP, CT and FT biomarkers in the calibrated population (blue histograms) and experimental ranges (black lines): TTF, time to peak force; FTD, force twitch duration; other abbreviations as defined in Table 1.

2.2 Ageing population of models and arrhythmogenic effects of age-related parameter modulation

The control model population was used to evaluate electro-mechanical response to agerelated modulations in four model parameters as reported in the literature: a reduction in the density of the potassium transient outward current Gto, maximal velocity of SERCA (Vmax_up), and an increase in the density of NCX current (K_NaCa) and CaL current (G_CaL). The reason for the parameter choice was experimental evidence on the age-related remodeling in these ionic mechanisms in different animal species [16-19] suggesting their possible role in human ageing myocardium. The direction of age-related parameter change was based on the available experimental data suggesting a reduction in the Gto and Vmax_up, and an increase in the K_NaCa and G_CaL [2].

Two protocols of age-related change in the model parameters were performed. First, we gradually changed each parameter individually by 20, 50, 70% in every model of the control population to evaluate the sensitivity of the entire population to the parameter variation and to compare the contribution of the parameters to the age-related abnormalities.

Then, we randomly generated 60 age-related sets of the 4 parameters using normal distributions for parameter deviation from the reference value with 25% mean and 10% SD and applied each set of the ageing parameters to every model in the control population. So we analyzed 60 aged sub-populations each consisted of 240 models and in total around 15’000 age-related models were evaluated in the trial.

Similarly, with virtual trials on arrhythmogenic drug-testing [4], we checked the occurrence of the age-related effects on the AP generation and Ca²⁺ handling associated with adverse events in the intact heart, particularly with the risk of TdP [14, 15].

We checked for repolarization and depolarization abnormalities (RA and DA respectively), as in [3, 4]. For models not displaying AP abnormalities, a more than 6% increase in APD₉₀: ΔAPD₉₀=(APD₉₀ aged-APD₉₀)/APD₉₀>6% was considered as risky for adverse consequences (adverse events). Moreover, we checked the implausible (ImP) model outputs out of the physiological ranges using the criteria in the lower part of Table 1. Such ImP outputs were usually associated with an increase in CT amplitude over an acceptable limit (see Table 1). The latter abnormalities always caused abnormalities in force/shortening generation. Models that do not fall into any of the mentioned above adverse categories, were labeled as “Norm”.

For the study, every aging model was paced at 1 Hz for 200 beats, for each parameter modulation. The last trace of each simulation was compared with the corresponding control (with no age-related parameter modulation) and automatically checked for the adverse events (ΔAPD>6%, RA, ImP).

3 Results

3.1 A human in silico control population of electro-mechanical models

Figure 1 demonstrates a summary for the behaviour of the control population of electromechanical models of human ventricular cardiomyocytes (Figure 1.A). In Figure 1. C, histograms of distributions of AP, CT and FT biomarkers versus the corresponding human experimental ranges are shown (see also Table 1). All the models show a normal AP phenotype, with all the AP and CT biomarkers within the experimental ranges.

The box plots for distributions of the model parameters in the control population (Figure 1. B) display a reduction in the non-implausible parameter range for V_{max_up}, G_CaL and P_NaK as compared with initial variation ranges due to filtering of implausible model samples. Moreover, it can be seen that both the medians for some parameters and the model traces significantly diverge from the reference values. Most of the models in the control population demonstrate higher membrane conductance for I_Na and I_CaL ionic currents, NCX and NKX exchangers, but lower conductance for I_Kr and reduced V_{max_up} as compared to the reference values.

3.2 Effects of individual modulation of age-dependent parameters

Here, we tested effects of age-dependent increase in G_CaL and K_NaCa, and a decrease in G_to and V_{max_up} as compared to the control parameter values in the population.

It is clearly seen in Figure 2, that the overall frequency of the adverse events (ΔAPD>6% + RA + ImP) increases with an increase in the age-dependent parameter deviation from the control for each tested parameter. At the same time, the ratio of individual risk indicators in the overall frequency of the adverse events is different between the model parameters.

The frequency of ImP models in the ageing population increases with an increase in the age-dependent parameter deviation for every parameter, where G_CaL increase has the highest impact on ImP occurrence in the model population. In the non-implausible ageing models which do not exhibit RA events, we showed an increase in both an average APD and variability in APD with an increase in age-related parameter deviation from the control value for each parameter varied (fig. 3). We found similar sensitivity of APD to individual agerelated change in G_CaL, K_NaCa and V_max against the control, while APD is almost insensitive to G_to, which mostly affects the phase of early repolarization of AP.

The fraction of Risky models with ΔAPD90>6% events notably increases with an agerelated increase in K_NaCa or a decrease in V_max against the control (fig. 2). An increase in K_NaCa causes the most notable increase in the fraction of ΔAPD₉₀>6% models. In opposite, the fraction of Risky models in the total ageing population decreases with an increase of GCaL over 50%.

The highest frequency of RA events in ageing models is shown in response to GCaL increase (fig. 2). While a number of RA models does not change essentially at increasing GCaL, their fraction in the non-ImP parameter space increases with parameter deviation from the control. As could be expected, Gto decrease did not affect APD, not induce RA and less frequently provokes ImP models as compared to other parameters varied. It is notable, that a 20% decrease in Gto caused a transition of approximately 20% of ageing models into the ImP category, whilst further decreasing Gto does not change the ImP fraction. We found that around 20% of models in the control population are near the “healthy” border in the parameter space, so even a small variation in a single parameter affects dramatically behaviour of such models.

3.3 Effects of combined modulation of age-dependent parameters

In the next series of simulations we randomly varied all the four age-dependent parameters simultaneously as a random 4-dimensional vector (see Methods). 60 sampled parameter vectors were generated and corresponding 60 different ageing populations based on the control population of models were computed. The random age-dependent parameter vector deviation D from the control (the distance from modified parameter vector to control one in the 4-dimensional parameter space) varied from 36% to 77% with average D of 54 ± 9%.

Figure 2 (bottom diagram) shows model classification in the 60 ageing populations at increasing distance D. Similarly with individual age-dependent parameter modulation, 4parameter vector deviations from the control throw out from 20% to 50% of ageing models from the non-implausible space, and the fraction of ImP models increases with increasing parameter deviation D (Fig. 2, right panels).

In the non-implausible ageing models with a combination of modified parameters, we also showed an increase in both an average APD and variability in ΔAPD with an increase in 4-parameter distance D from the control (fig. 3). We found positive correlations between the ageing parameter deviation D and the frequency of the RA events in the ageing population of models. At the same time, at the combined parameter modulation the frequency of ageing models with ΔAPD90>6% decreases within the non-implausible parameter space (fig. 2).

Fig. 2. A Classification of models adverse events during aging. Colour maps indicate “Norm” (green), “Risky” (blue tones), “RA” (orange), and “ImP” (rose) classes in the ageing populations. “Risky” class contains models showing ΔAPD90>6% events. Top frequency diagrams show results of age-related deviation of individual parameters, bottom diagram shows classification for 60 ageing populations with 4 parameter varied. Dashed line shows smoothed borderline between classes. B Fraction of ImP models and arrhythmogenic score dependencies on the age-related parameter deviation. Orange, green, red and purple lines show dependencies at GCaL, KNaCa, Gto and Vmax_up individual modulation, respectively. Black dots show scatter plots for 60 ageing models with blue linear regression lines (p < 0.01).

Fig. 3. Effects of modulation of age-dependent parameters on APD. Blue and orange lines show APD90 and APD50 values respectively, plus or minus standard deviations are shown as semitransparent overlays. Left figure shows results of 60 ageing populations with simultaneously varied parameters, right figures show results of age-related deviation of individual parameters.

4 Discussion and conclusions

In this study, we developed a new simulation-based approach to predict the consequences of age-related cellular remodeling. Currently, a large experimental data on myocardial remodeling associated with ageing is known, including ionic remodeling affecting the function of cardiomyocytes [1, 2]. In this study, we focused on in-silico assessment of electrophysiological consequences of the age-related modulations in the cellular ionic parameters and on the prediction of the arrhythmogenic risk of the cellular remodeling in aging using human ventricular cellular models of electro-mechanical activity. First, we developed a control population of electro-mechanical models. We used the same acceptable biomarker ranges for model calibration as proposed in [4] with additional restrictions on the CT amplitude and the mechanical function of the cell based on the Frank-Starling law [12]. Moreover, we extended acceptable range for CTD₅₀ based on data from [11] (see Table 1). With using those restrictions, the model acceptance was rather low with about 240 from 20000 calculated models met all filters for AP, CT and FT biomarkers. Such a low acceptance level reflects an essential role of the CT and mechanical biomarkers in the model selection.

To simulate ageing populations, we randomly selected 60 sets of age-related parameter modification by uni-directional deviation of the 4 ionic parameters from the control values:

G_CaL and K_NaCa were gradually increased in aging, while V_{max_up} and G_to were reduced according to the available experimental data [16-19]. These parameter modulations were applied to each of the models from the control population, and age-related changes in the model outputs were evaluated in the populations. It is still difficult to evaluate quantitatively how large the molecular changes in age-related conductivity are in human cardiomyocytes. That is why we simulated parameter deviations in a rather wide range, to be able to reveal possible age-dependence of the effects.

Our population of models was able to reproduce aging effects on cellular biomarkers observed in mammals. That is a prolongation of AP and CT [1, 2], particularly reflected in an increase in QT intervals in ECG. In concordance with experimental findings, our population of ageing models showed an increase in mean APD (fig. 3) and CT duration (not shown) with increasing age-related parameter deviation from the control.

We used frequency of RA events and an increase in APD over threshold level (ΔAPD>6% as suggested in [3, 4]) as possible scores of arrhythmogenicity of the agedependent ionic remodeling in human cardiomyocytes. The population of ageing models predicted an increase in the frequency of the RA events at both individual and combined parameter deviations. For ΔAPD>6% score the results are less straightforward. This index increases with individual age-dependent parameter deviation from control except of large deviations of GCaL where the frequency of RA events becomes higher. At the same time, for combined parameter deviation the negative regression was found between the frequency of ΔAPD90>6% models in non-implausible ageing population and parameter distance D. The mechanisms of this inconsistency have to be further analysed.

The population based analysis allowed us to compare the contribution of different excitation-contraction mechanisms in age-related effects on cellular biomarkers. An increase in the ICaL current, Na⁺-Ca²⁺ exchanger activity and a decrease in the Ca²⁺ SERCA uptake are predicted to be significant factors of age-related cellular remodeling increasing risk of arrhythmogenic abnormalities in elderly. The essential contribution of Vmax reduction in the tested adverse electrical events in ageing population is not obvious because this parameter does not affect the membrane conductance directly. It affects CT time course and duration which causes modulation of Ca²⁺ dependent currents resulting in APD prolongation over the threshold and RA occurrence.

In conclusion, our main findings are as follows:

1. We developed a control population of electro-mechanical cellular models of humanventricular cardiomyocytes which can be used for in-silico studies of the effects of different cellular interventions, including drug-testing, age-related remodeling, cellular pathological abnormalities.

2. The use of arrhythmogenic scores based on the RA and ΔAPD > 6% criteria allow one to predict and classify pro-arrhythmic risk of age-dependent remodeling of intracellular mechanisms of excitation-contraction coupling.

This work supported by RBRF # 18-015-00368, IIF UrB RAS theme #AAAA- A18118020590031-8, RF Government Act #211 of March 16, 2013.

References

All Tables

All Figures

Fig. 1. A control population of human ventricular electro-mechanical models. A Traces of action potential (AP), Ca2+ transient (CT) and force twitch (FT) in the 240 models in the calibrated population. B Ionic parameter profile of the calibrated population, shown as scaling factors with respect to the baseline model. Central mark at box-plots corresponds to the median, box limits are the 25th and 75th percentiles, and whiskers extend to the most extreme data points. C Distributions of AP, CT and FT biomarkers in the calibrated population (blue histograms) and experimental ranges (black lines): TTF, time to peak force; FTD, force twitch duration; other abbreviations as defined in Table 1.

In the text

Fig. 2. A Classification of models adverse events during aging. Colour maps indicate “Norm” (green), “Risky” (blue tones), “RA” (orange), and “ImP” (rose) classes in the ageing populations. “Risky” class contains models showing ΔAPD90>6% events. Top frequency diagrams show results of age-related deviation of individual parameters, bottom diagram shows classification for 60 ageing populations with 4 parameter varied. Dashed line shows smoothed borderline between classes. B Fraction of ImP models and arrhythmogenic score dependencies on the age-related parameter deviation. Orange, green, red and purple lines show dependencies at GCaL, KNaCa, Gto and Vmax_up individual modulation, respectively. Black dots show scatter plots for 60 ageing models with blue linear regression lines (p < 0.01).

In the text

	Fig. 3. Effects of modulation of age-dependent parameters on APD. Blue and orange lines show APD90 and APD50 values respectively, plus or minus standard deviations are shown as semitransparent overlays. Left figure shows results of 60 ageing populations with simultaneously varied parameters, right figures show results of age-related deviation of individual parameters.
In the text