A Novel In Silico Electromechanical Model of Human Ventricular Cardiomyocyte

Contractility has become one of the main readouts in computational and experimental studies on cardiomyocytes. Following this trend, we propose a novel mathematical model of human ventricular cardiomyocytes electromechanics, BPSLand, by coupling a recent human contractile element to the BPS2020 model of electrophysiology. BPSLand is the result of a hybrid optimization process and it reproduces all the electrophysiology experimental indices captured by its predecessor BPS2020, simultaneously enabling the simulation of realistic human active tension and its potential abnormalities. The transmural heterogeneity in both electrophysiology and contractility departments was simulated consistent with previous computational and in vitro studies. Furthermore, our model could capture delayed afterdepolarizations (DADs), early afterdepolarizations (EADs), and contraction abnormalities in terms of aftercontractions triggered by either drug action or special pacing modes. Finally, we further validated the mechanical results of the model against previous experimental and in silico studies, e.g., the contractility dependence on pacing rate. Adding a new level of applicability to the normative models of human cardiomyocytes, BPSLand represents a robust, fully-human in silico model with promising capabilities for translational cardiology.

1 Introduction

The future of diagnosis and treatment in cardiology progressively depends on advanced methods in imaging, gene profiling, and pharmaceutical technologies. Despite the recent advances in health technologies, the current empirical clinical investigations face serious challenges as the complexity of therapeutic interventions, prognosis, and the possibility of classifying the treatments grow. Specifically, identifying the optimal treatment strategy with a degree of statistical significance poses serious challenges to current empirical routes in cardiology (Niederer et al., 2019). Furthermore, as precision medicine emerges (Forouzandehmehr et al., 2022), the proven pathophysiological variability between individuals highly augments the detail in the diagnostic process and data, thus, finding an optimal patient-specific solution becomes increasingly difficult (Niederer et al., 2019). Cardiac computational models, based on established physiological and engineering principles, offer a capable framework that not only can be fed by large datasets but also enable mechanistic and integrative simulations leading to disclose novel insights in cardiac pathophysiology (Niederer et al., 2019).

Early physiologically constrained computational models of cardiac cells could quantitatively translate the protein functions into developing cellular phenotypes (Niederer et al., 2019). During the past decade, these models have also incorporated functional characteristics of ion channels, cellular pumps, transporters, and buffers making them promising candidates in preclinical studies (O’Hara et al., 2011; Tomek et al., 2019; Bartolucci et al., 2020; Paci et al., 2021). Currently, as the cardiac contractility data become increasingly available, together with new recording techniques (Ahola et al., 2018), mathematical models of cardiomyocytes (CMS) are developed to predict dynamics of contraction combined with simulation of drug effects alongside the electrophysiology. Toward building models of myocyte electromechanics, elements of contractility have been developed with different levels of complexity focusing on animals (Rice et al., 2008; Campbell et al., 2010; Sheikh et al., 2012; Land et al., 2013) and human CMs (Land et al., 2017). Initial efforts on developing reliable models to capture the electromechanics of human adult CMs have been initiated recently (Lyon et al., 2020; Margara et al., 2021). Margara et al., integrated an established human-based developed contractile element (Land et al., 2017) into the gold standard in silico model of human ventricular cell electrophysiology (O’Hara et al., 2011) and into their new ToR-ORd model (Tomek et al., 2019) to predict ventricular active tension generation alongside action potential (AP) and calcium transients (CaT). Also Lyon et al. complemented the ORd model with a contractile element: their choice was the MedChem model of sarcomere mechanics (Dupuis et al., 2018), which they used to assess the impact of β-adrenergic stimulation and sarcomere length on CaTs and force (Lyon et al., 2020).

Our recently published BPS2020 model of the human adult ventricular CM electrophysiology (Bartolucci et al., 2020), holds significant improvements compared with the original ORd model (e.g., the simulation of the experiments with the correct extracellular K⁺ concentration used in vitro or the generation of DADs) and includes specific mechanisms not simulated by ToR-ORd (e.g., the inverse APD₉₀-[Ca²⁺]_o dependence). Given these improvements in simulating electrophysiology phenomena, it is worth investigating how we can expand the spectrum of BPS2020 simulation, by using it as the basis for a new electromechanical human CM in silico model.

We have integrated one of the most recent human contractile machinery (LandCE) (Land et al., 2017) into BPS2020 (Bartolucci et al., 2020). As done in Margara et al. (2021), we chose LandCE as it is a contractile element validated against human data. Our goal was to investigate the capabilities of this newly integrated electromechanical model, BPSLand, by evaluating active tension generation and contractility abnormalities e.g., aftercontractions, that can be activated either by drug action or special pacing conditions. BPSLand is a robust, fully-human, in silico model meeting the computational expectations from both departments, the electrophysiology and contractility, with the potential for facilitating the translation of biophysical and pharmacological functions into pre-clinical readouts.

2 Methods

2.1 In Vitro Data

To calibrate the BPSLand model, we used a dataset of isometric active tension (Ta) biomarkers recorded from human isolated ventricular CMs (Land et al., 2017; Margara et al., 2021) and a dataset of action potential (AP) biomarkers from human isolated ventricular endocardial CMs (O’Hara et al., 2011; Bartolucci et al., 2020). The Ta biomarker dataset includes measures from strips of the left ventricular myocardium (Mulieri et al., 1992), left ventricular trabeculae (Pieske et al., 1996) and right ventricular trabeculae (Rossman et al., 2004) (additionally considered by Margara et al., 2021). Both datasets were recorded at 37°C. Ta biomarkers are the peak tension (TaPeak), the relaxation time at 50% and 95% (TaRT₅₀, TaRT₉₅) and the time-to-peak (TaTTP). AP biomarkers are the duration at 30%, 40%, 50%, 70% and 90% (APD_30, APD₄₀, APD₅₀, APD₇₀, APD₉₀), the maximum upstroke velocity (dV/dt_max), the peak voltage (VPeak) and the resting membrane potential (RMP). In silico biomarkers were computed as in Margara et al. (2021). As we previously reported in Bartolucci et al. (2020), we simulated the AP biomarkers for calibration at [K⁺]_o = 4 mM. Conversely, as no information was reported on the in vitro Ta biomarker ranges, we run our simulations at the standard concentration [K⁺]_o = 5.4 mM.

To validate BPSLand, we used the following human data: 1) APD rate dependence, restitution and accommodation data in control condition and with current blockers from endocardial CMs (O’Hara et al., 2011) (see Supplementary Methods for details and Supplementary Table S1); 2) TaPeak, TaRT₅₀ and CaT relaxation time at 50% (CaRT₅₀) rate adaptation data (Pieske et al., 1995; Janssen and Periasamy, 2007); 3) TaPeak transmural heterogeneity data from sub-epicardial, mid-myocardial, and sub-endocardial specimens (Haynes et al., 2014).

2.2 Integration of the Land Contractile Element Into the BPS2020 Model

The original BPS2020 model (Bartolucci et al., 2020) was based on the seminal O'Hara-Rudy model of the human ventricular AP (O’Hara et al., 2011) and it features two cytosolic compartments, the subspace and the bulk myoplasm, and the sarcoplasmic reticulum (SR) represented with a single compartment. It includes the following ion currents: fast and late Na⁺ currents (I_NaF, I_NaL), transient outward K⁺ current (I_to), L-type Ca²⁺ current (I_CaL), also with its Na⁺ and K⁺ components (I_CaNa, I_CaK), the rapid, slow and inward rectifying K⁺ currents (I_Kr, I_Ks, I_K1), the Na⁺/Ca²⁺ exchanger divided in its cytosolic and subspace components (I_NCXi, I_NCXss), the Na⁺/K⁺ pump (I_NaK), Na⁺, K⁺ and Ca²⁺ background currents (I_Nab, I_Kb, I_Cab) and the sarcolemmal Ca²⁺ pump (I_pCa). Ca²⁺ fluxes from/to SR are the RyR-sensitive Ca²⁺ release (J_rel), the SERCA pump (J_up) and a leakage flux (J_leak).

We integrated LandCE into BPS2020 following the approach presented in (Margara et al., 2021). Shortly, LandCE takes as input the intracellular Ca²⁺ concentration [Ca²⁺]_i computed by BPS2020, to update a new state variable CaTRPN, representing the fraction of troponin C units which bound to Ca²⁺.

$\frac{dCaTRPN}{dt}= k_{TRPN}\left({\left(\frac{{\left[C{a}^{2+}\right]}_i}{{\left[C{a}^{2+}\right]}_{T\text{50}}}\right)}^{n_{TRPN}}\left(\text{1}-CaTRPN\right)-CaTRPN\right)$

BPS2020 receives as feedback the amount of Ca²⁺ buffered by troponin C, [Ca²⁺]_TRPN, to update the intracellular Ca²⁺ concentration.

$\frac{d{\left[C{a}^{\text{2}+}\right]}_i}{dt}= {\mathbf{\beta }}_{Cai}\left(-\left(I_{pCa}+I_{Cab}-\text{2}{I}_{NaCa,i}\right)\frac{A_{cap}}{\text{2}F{v}_{myo}}-J_{up} \frac{v_{sr}}{v_{myo}}+ J_{diff,Ca} \frac{v_{ss}}{v_{myo}}- \frac{d{\left[\mathbf{E}\mathbf{G}\mathbf{T}\mathbf{A}\right]}_i}{dt}- \frac{d{\left[C{a}^{\text{2}+}\right]}_{TRPN}}{dt}\right)$

$\frac{d{\left[C{a}^{\text{2}+}\right]}_{TRPN}}{dt}={\left[C{a}^{\text{2}+}\right]}_{TRPN,max}\frac{dCaTRPN}{dt}$

where ${\left[C{a}^{2+}\right]}_{TRPN,max}$ represents the maximum Ca²⁺ concentration that can bind to troponin C.

2.3 Optimization of the BPSLand Model

The structure of the cost function used for both optimizations is the same as in Paci et al. (2018b)

$Cost={\displaystyle \sum _{\mathbf{1}}^{N_{biom}}{w}_i\ast Cos{t}_i}$

$\mathit{Cos}{\mathit{t}}_{\mathit{i}}=\frac{\left({\mathit{b}}_{\mathit{i},\mathit{sim}}<\mathit{ L}{\mathit{B}}_{\mathit{i}}\right){\left({\mathit{b}}_{\mathit{i},\mathit{sim}}-\mathit{L}{\mathit{B}}_{\mathit{i}}\right)}^2+\left({\mathit{b}}_{\mathit{i},\mathit{sim}}>\mathit{ U}{\mathit{B}}_{\mathit{i}}\right){\left({\mathit{b}}_{\mathit{i},\mathit{sim}}-\mathit{U}{\mathit{B}}_{\mathit{i}}\right)}^2}{\mathbf{0.5}(\mathit{L}{\mathit{B}}_{\mathit{i}}+\mathit{U}{\mathit{B}}_{\mathit{i}})}$

where b_i,sim is the ith simulated biomarker, LB_i the ith experimental lower bound for b_i,sim, UB_i the ith experimental upper bound for b_i,sim, w_i the weight for each biomarker’s cost (Supplementary Table S2) and N_biom the number of biomarkers used for optimization. Briefly, if the simulated ith biomarker is smaller than LB_i or greater than UB_i, the error is computed as the squared distance between the simulated biomarker and the bound, normalized by the center of mass of [LB_i, UB_i]. Finally, in order to minimize the active tension T_a, we included one additional term to Cost, obtaining the final cost function

${\mathit{Cost}}_{\mathit{TOT}}={\mathit{w}}_{\mathit{minTa}}\mathit{\ast }\mathbf{min}\left({\mathit{T}}_{\mathit{a}}\right)+\mathit{Cost}$

with $w_{minTa}$ the weight of the minimun active tension.

2.3.1 Step 1: Hybrid Optimization on the LandCE Parameters

After integrating LandCE into BPS2020, we first optimized the LandCE parameters using a hybrid approach combining first a genetic optimization (Matlab function ga), followed by the simplex optimization [Matlab function fminsearchbnd (D’Errico, 2022)]. The parameters optimized in this first step are only the LandCE parameters listed in Table 1. The optimization ranges for the LandCE parameters are the same as in the original LandCE publication (Land et al., 2017), except for the tropomyosin Ca²⁺ sensitivity ([Ca²⁺]_T50), for which we chose [0.5, 0.6] instead of [0.8, 0.9]. As the original range [0.8, 0.9] increased substantially the CaT peak, i.e. the systolic Ca²⁺ (+22%), we decided not to affect the BPS2020 electrophysiology and chose [0.5, 0.6] as it preserved the original BPS2020 CaT peak.

TABLE 1

TABLE 1. Contractility and electrophysiology biomarkers used for the BPSLand optimization, with their ranges.

For this first optimization step, we considered five contractility and two electrophysiology biomarkers: active tension peak (TaPeak), time-to-peak (TaTTP), relaxation time to 50% and 95% of the diastolic level (TaRT₅₀ and TaRT₉₅) and the minimum of the diastolic active tension, systolic and diastolic intracellular free Ca²⁺ (CaSys and CaDiast). The acceptable ranges for these biomarkers were taken from the original Land publication (Land et al., 2017) for TaTTP, TaRT₅₀ and TaRT₉₅, from Margara et al. (2021) for TaPeak, while we chose to set the ranges for CaSys and CaDias as ±5% of their original values (Bartolucci et al., 2020), in order to keep the electrophysiology the most similar to the original BPS2020 model. At the end of this first step, we obtained an electromechanical model whose electrophysiology biomarkers were not significantly affected by the LandCE and correctly simulated TaRT₉₅ while the remaining contractility biomarkers were close to their respective lower bounds.

2.3.2 Step 2: Simplex Only

In order to capture the remaining contractility biomarkers, we then run a second simplex optimization on the Ca²⁺ fluxes of the SERCA pump (J_up) and the RyR-sensitive release (J_rel), using all the constraints in Table 2, and additional constraints on the AP biomarkers. In particular, for resting potential (RMP), peak voltage (VPeak), maximum upstroke velocity (dV/dt_max), AP duration at 40%, 50% and 90% (APD₄₀, APD₅₀ and APD₉₀), and the triangulation metric (Tri₉₀₄₀), we set the lower and upper bounds as ±5% of their values in the original BPS2020 model, which were fit against the experimental data (O’Hara et al., 2011) at [K⁺]_o = 4 mM. We chose these parameters as we did not want to change the ion current parameters of BPS2020, derived from the ORd model and partially fit experimental data in Bartolucci et al. (2020). As ranged for manually tune J_rel,max and J_up,max, we chose ±20% of their original values 20e-3 (1/ms) and 3.13 (mM/ms), respectively. At the end of this second step, the electromechanical model correctly simulated TaRT₉₅ and TaPeak while the remaining contractility biomarkers were close to their respective lower bounds. However, we missed one of the key features of BPS2020, i.e., the inverse relationship [Ca²⁺]_o – APD₉₀, which was otherwise simulated at the end of the first step.

TABLE 2

TABLE 2. In vitro contractility and electrophysiology biomarkers used in the cost function and their goal ranges.

2.3.3 Step 3: Manual Tuning

In order to restore the [Ca²⁺]_o – APD₉₀ relationship, we added one final step to our pipeline, where we did a minor manual re-tuning of J_rel,max (0.0240 →0.0220 1/ms) and J_up,max (3.1333 → 3 mM/ms), still considering their lower and upper bounds as in Table 1. The final model is named BPSLand and its parameter values are reported in Table 3. Supplementary Tables S2–S4 summarize the weights, parameters and biomarkers obtained after each of the three optimizations steps. Supplementary Table S5 shows the impact of the manual tuning of J_rel,max and J_up,max on the [Ca²⁺]_o – APD₉₀ relationship.

TABLE 3

TABLE 3. Final BPSLand parameter set.

2.3.4 Rate Dependence

To test the active tension dependence on the applied pacing rate, we paced BPSLand at 0.5, 1, 1.5, 2, 2.5 and 3 Hz for 1,000 beats to reach the steady state, using [K⁺]_o = 5.4 mM, [Ca²⁺]_o = 1.8 mM and [Na⁺]_o = 144 mM as extracellular ion concentration. We then compared simulated TaPeak, TaRT₅₀ and CaRT₅₀ with the in vitro data by Pieske et al., (1995) and Janssen and Periasamy (2007).

2.3.5 Heterogeneity

To simulate transmural heterogeneity, i.e., simulating epicardial (EPI) and mid-myocardial (M) CMs in addition to endocardial (ENDO), we used the same scaling factors reported in Bartolucci et al. (2020) for I_NaL, I_to, I_CaL, I_Kr, I_Ks, I_K1, I_NCX, I_NaK, I_Kb, J_rel and J_up (Supplementary Table S6).

2.3.6 EAD and DAD Simulations to Trigger Active Tension Abnormalities

We assessed the occurrence of early-afterdepolarization (EADs) and aftercontractions in the BPSLand following three different protocols. First, we simulated the administration of quinidine considering the drug effects on I_Na, I_Kr, I_CaL, I_Ks and I_to, using the IC₅₀ and Hill’s coefficients reported in Passini et al. (2017) and Paci et al. (2018a) and the single pore block model, as in Paci et al. (2021) (Supplementary Table S7). We tested three drug concentrations, namely 10, 15 and 20 μM at the standard extracellular ion concentrations ([K⁺]_o = 5.4 mM, [Ca²⁺]_o = 1.8 mM, [Na⁺]_o = 144 mM) and cycle length (CL) of 4,000 ms. The second EAD protocol simulated dofetilide, similarly to what we did in Bartolucci et al. (2020). Shortly, we simulated the administration of 0.1 µM dofetilide at CL = 5,000 ms and extracellular concentrations experimentally used by Guo et al. (2011) ([K⁺]_o = 5 mM, [Ca²⁺]_o = 2 mM, [Na⁺]_o = 137 mM), using the I_Kr drug binding values reported by Dutta et al. (2017). We simulated quinidine and dofetilide effects on the endocardial BPSLand model and we anticipate no EADs nor aftercontractions, despite the remarkable AP prolongation. Conversely, the same tests performed on the M cell version, resulted in EADs and aftercontraction.

To assess the occurrence of delayed afterdepolarizations (DADs) we used the same protocol as in Li and Rudy (2011): we fast paced BPSLand for 1,500 beats (BCL = 275 ms) and then we triggered one long beat (BCL = 10,000 ms).

3 Results

3.1 The BPSLand Model

We report the AP, [Ca²⁺]_i, [Ca²⁺]_ss and Ta traces simulated at [K⁺]_o = 5.4 mM in Figure 1, together with the simulations for variable [Ca²⁺]_o to highlight the inverse APD₉₀-[Ca²⁺]_o dependence, which was described first by Severi et al. (2009) and then observed in vitro and in vivo (Leitch, 1996; Nagy et al., 2013), but failed to be simulated by many in silico models, including the original ORd (O’Hara et al., 2011) and ToR-ORd (Tomek et al., 2019). In details, for increasing [Ca²⁺]_o = 0.9, 1.8 and 2.4 mM, APD₉₀ equals to 251.4, 239.9 and 237.1 ms. Table 4 reports the AP and Ta biomarkers, with the in vitro ranges used for the BPSLand calibration, together with additional CaT biomarkers: CaT duration at 50% and 90% (CTD₅₀, CTD₉₀) and amplitude (CaAmp). All the AP biomarkers are within the experimental ranges, as well as TaPeak and TaRT₉₅. Conversely, TaRT₅₀ and TaTTP are very close to their respective experimental lower bounds, although out of the in vitro ranges. The comparison of AP, Ca²⁺ and Ta traces simulated with [K⁺]_o = 4 mM and [K⁺]_o = 5.4 mM is presented in Supplementary Figure S1.

FIGURE 1

FIGURE 1. Illustrative traces simulated by BPSLand ([K⁺]_o = 5.4 mM). (A) Action potential. (B) Cytosolic Ca²⁺ concentration. (C) Subspace Ca²⁺ concentration. (D) Active tension. (E) Inverse action potential duration dependence on the extracellular Ca²⁺ concentration.

TABLE 4

TABLE 4. The electrophysiology and contractility biomarkers simulated by the original BPS2020 and the new BPSLand models, compared to in vitro data.

3.2 Electrophysiology and Contractility Dependence on Pacing Rate

The APD rate adaptation tests reported in Bartolucci et al. (2020) were repeated using BPSLand, to show that introducing LandCE did not affect the capability of the new model in simulating the old data. Briefly, BPSLand simulated the in vitro data as satisfactorily as BPS 2020, outperforming the original ORd model (Supplementary Figures S2, S3).

In this section we validate BPSLand against two additional in vitro datasets of rate adaptation of TaPeak, TaRT₅₀ and CaRT₅₀, not considered for BPS2020. Figure 2A shows the qualitative agreement of our model with the data published by Janssen and Periasamy (2007) in terms of TaPeak-pacing rate dependence. In particular, we successfully simulate the linearity of such dependence. In Figure 2B, we considered the rate dependence of TaRT₅₀, considering in vitro data by Pieske et al., 1995 and Janssen and Periasamy (2007). BPSLand simulations qualitatively agree both with the Janssen07 and the Pieske95 experiments, although TaRT₅₀ is lower at the slowest pacing rates. This discrepancy is due to the TaRT₅₀ in vitro range used to calibrate the BPSLand model at 1 Hz, i.e. [109, 125] ms (purple line). BPSLand is positioned at the interval lower bound (108.4 ms), while Janssen07 data at the upper bound (125 ms) and Pieske95 is out of bound (137.2 ms). Conversely, BPSLand shows quantitative agreement with the Pieske95 CaRT₅₀ data (Figure 2C). A comparison with the ToR-ORd+Land model is shown in Supplementary Figure S4. Furthermore, the length dependence properties of the BPSLand model is presented in Supplementary Figure S5 and was performed in the same way proposed by Margara et al. (2021) in their original Supplementary Section S6.

FIGURE 2

FIGURE 2. Active tension dependence on pacing rate and comparison with the in vitro data from Pieske et al. (1995) and Janssen and Periasamy (2007). (A) Normalized active tension peak. (B) Active tension relaxation time. (C) Ca²⁺ transient relaxation time.

3.3 Transmural Heterogeneity

Figure 3 shows how BPSLand simulates the transmural heterogeneity in terms of electrophysiology and contractility. Our simulations are in agreement with the ToR-ORd+Land and ORd+Land models presented in Margara et al. (2021). In terms of APD, the M model has the longest APs, followed by the ENDO and EPI models. In terms of CaTs and active tension, the M model shows the highest peaks, followed by EPI and ENDO. Haynes et al. reported transmural heterogeneity data of isometric active tension peaks in human heart preparations, showing similar average active tension in EPI and ENDO preparations (although EPI < ENDO), and greater in M specimens (Haynes et al., 2014). We simulate an EPI TaPeak (17.4 kPa) slightly greater than ENDO (15.6 kPa), while the M model produces greater TaPeak (34.8 kPa). This is the same trend simulated by the ToR-ORd+Land model (TaPeak M > EPI > ENDO), although the absolute TaPeak values are considerably greater in ToR-ORd+Land than in BPSLand. As in (Margara et al., 2021) the authors suggested that the Ca²⁺ sensitivity in ENDO CMs could be higher than in EPI cells, we tested how much upscaling of [Ca²⁺]_T50 is required in EPI BPSLand to bring the simulated EPI TaPeak even closer to the experiments (Haynes et al., 2014). In fact, [Ca²⁺]_T50 is not considered as one of the parameters to change when switching from ENDO to EPI models. The purple star in Figure 3E show that a ×1.1 upscale produces an EPI TaPeak matching the experiments. A comparison of the transmural heterogeneity with the ToR-ORd+Land model is also reported in Supplementary Figure S6.

FIGURE 3

FIGURE 3. Transmural heterogeneity simulations with BPSLand in endocardial (green), epicardial (blue) and mid-myocardial (red) modes. (A) Action potentials. (B) Cytosolic Ca²⁺ concentration. (C) Active tension. (D–E) Comparison of the simulated and experimental active tension peak magnitude across endocardial, epicardial and mid-myocardial cell types. The simulated purple star represents an additional epicardial simulation where we tested a small increments of the calcium sensitivity, ×1.1 the baseline [Ca²⁺]_T50 value, to obtain an active tension even closer to the experimental data (Haynes et al., 2014).

3.4 EADs, DADs and Aftercontractions

The endocardial BPSLand model did not produce EADs just by administering quinidine or dofetilide, despite the extreme APD₉₀ prolongation up to +272% with 0.1 µM dofetilide; +398%, +489%, +563% with the three increasing quinidine doses.

Conversely, the M cell version, characterized by smaller I_Kr and larger I_CaL, reacted to both drugs with EADs and, in some cases, aftercontractions. The simulations shown in Figure 4 are noteworthy: for both the intermediate and high quinidine concentrations (15 µM in the second row and 20 µM in the third row) quinidine triggers EADs, but only some of them have a correspondent aftercontraction. This is due to the different mechanisms underlying each EAD and it is well summarized in case of 20 µM quinidine, reported in more detail in Figure 5. The smaller EADs due to I_CaL reactivation (I_CaL-driven), e.g., t ∼ 1.3 s or t ∼ 5.3 s do not have a corresponding aftercontractions. On the other hand, other EADs are triggered by a spontaneous Ca²⁺ release from the SR through J_rel, e.g., t ∼ 2.2 s or t ∼ 10.2 s, which pours into the cytosol enough Ca²⁺ to trigger the contractile element to produce an aftercontraction. Therefore, from this result, we can hypothesize there is not a 1:1 EAD-aftercontraction correspondence, since aftercontractions require enough Ca²⁺ to start, as in the case of J_rel intervention.

FIGURE 4

FIGURE 4. Illustrative traces of the membrane potential (left column) and active tension (right column) simulated by the M cell BPSLand with low (10 µM), intermediate (15 µM) and high (20 µM) quinidine concentrations. The intermediate and high doses trigger early afterdepolarizations and aftercontractions.

FIGURE 5

FIGURE 5. Early afterdepolarizations (EADs) triggered by 20 µM of quinidine and their underlying mechanisms in the M cell BPSLand. The smaller EADs due to I_CaL reactivation (I_CaL-driven), e.g., t ∼ 1.3 s or t ∼ 5.3 s do not have a corresponding aftercontraction. Conversely, EADs trigger by a spontaneous Ca²⁺ release from the sarcoplasmic reticulum through J_rel, e.g., t = ∼ 2.2 s or t = ∼ 10.2 s, have a corresponding aftercontraction, since J_rel pours into the cytosol enough Ca²⁺ to trigger the contractile element.

We observed a similar result with dofetilide in Figure 6 where the dofetilide simulation resulted in EADs and aftercontractions. Also in this case, the EADs are triggered by spontaneous Ca²⁺ release from the SR through J_rel, as shown in the third panel.

FIGURE 6

FIGURE 6. Illustrative early afterdepolarizations and aftercontractions triggered by 0.1 µM of dofetilide in the M cell BPSLand.

Following the DAD Li et al. protocol, BPSLand triggered an unpaced beat further followed by several DADs (Figure 7). The fast pacing protocol led to the accumulation of Ca²⁺ in the SR (oscillations in [1.76, 2.12] mM instead of [1.20, 1.47] mM), which was spontaneously released by J_rel during the diastolic phase of the last long beat. These unpaced releases of sarcoplasmic Ca²⁺ not only triggered the anticipated AP and DADs (as we already showed in Bartolucci et al., 2020), but it also was enough to trigger aftercontractions (Figure 7D inset).

FIGURE 7

FIGURE 7. Aftercontractions triggered by anticipated beats and DADs in the endocardial (left column) and M cell (right column) BPS 2020. The action potential at t = 0 is the long beat at BCL = 10,000 ms, following 1,500 beats at BCL = 275 ms. The action potential at t ∼ 6 s (left) and t ∼ 4.2 s/4.7 s/6.4 s (right) is triggered by the spontaneous Ca²⁺ release from the sarcoplasmic reticulum and not by external pacing. (A) Membrane potential. (B) Ca²⁺ release flux from the sarcoplasmic reticulum. (C) Cytosolic Ca²⁺ concentration. (D) Active tension with aftercontractions. The zoomed inset on the left column highlights the small aftercontractions corresponding to the DADs following the anticipated action potential.

4 Discussion

In this work, we present an updated version of our BPS2020 model of the human ventricular AP (Bartolucci et al., 2020), that we enhanced with the contractility model presented by Land et al. (2017). The potential of in silico models is getting more and more recognition both by industry and regulators for specific applications, e.g., cardiac safety pharmacology (Li et al., 2020; Musuamba et al., 2021). However, most of the current cardiac cell models focus mainly on electrophysiology, i.e., AP and Ca²⁺ handling, not considering the fact that the heart behaves like a pump, and therefore the contractile activity of CMs is surely worth of interest. Most of the diseases of interest modelled so far within silico CM models mainly affected specific ion channels [long QT syndrome (Clancy and Rudy, 2002; Paci et al., 2017, 2018a; Kernik et al., 2020)] or Ca²⁺ handling [catecholaminergic polymorphic ventricular tachycardia (CPVT) (Koivumäki et al., 2018)]. Conversely, hypertrophic cardiomyopathy (HCM), the most widespread genetic cardiac disorder, primarily associates with pathogenic variants in protein genes of sarcomere (Santini et al., 2020). In fact, most of pathogenic variants in HCM are hosted by myosin binding protein C and adult cardiac myosin isoforms that are mainly programmed by MYBPS3 and MYH7 genes, respectively, (Toepfer et al., 2020). These variants are responsible for myocardium hypercontractility (Sarkar et al., 2020), impaired contractile relaxation (Toepfer et al., 2020), arrhythmogenesis, diastolic dysfunction and heart failure (Sarkar et al., 2020). Furthermore, the hypoxia-induced lack of oxygenation in ischemia impairs the orchestrate of molecular events leading to normal ventricular contraction (Katz, 1973). Finally, the glycation of myofilaments in diabetes, a major risk factor in heart failure, correlates with significant reduction in calcium sensitivity of the sarcomere (Papadaki et al., 2022) that cannot be captured in electrophysiology-only models. The same applies to new drugs directly targeting sarcomere dynamics, e.g., blebbistatin, omecamtiv mecarbil and mavacamten (Rahman et al., 2018; Awinda et al., 2020; Fülöp et al., 2021): with no in silico contractility description, it is not possible to properly simulate them.

Therefore, the goal of our work was to provide the new and validated BPSLand in silico model of human adult CMs, combining both electrophysiology and contractility. As the electrophysiology description by BPS2020 carried a few novelties, especially the APD-[Ca²⁺]_o relationship and an extended and more reliable description of Ca²⁺ handling, including the generation of DADs, it was important to us to create a model able to translate such novelties also to contractility. Of note, we did not aim to simulate here specific pathological conditions affecting contractility, as this will be topic for future works.

4.1 Development of the BPSLand Model

We followed the same strategy published by Margara et al. (2021) for their ToR-ORd+Land model, to integrate the electrophysiology described by BPS2020 and the contractility of LandCE: as forward mechanism, LandCE takes as input the cytosolic Ca²⁺ concentration computed by BPS2020, to compute the fraction of troponin C units bound to Ca²⁺, and this new flux of Ca²⁺ towards the sarcomere is then included in the equation regulating the BPS2020 cytosolic Ca²⁺ concentration, to close the loop. In terms of mathematical formulation, the process was straightforward, as BPS2020 and ToR-ORd are both based on the original ORd model. For the optimization of the model, we built our cost function with the same biomarkers (TaPeak, TaTTP, TaRT₅₀ and TaRT₉₅) and experimental ranges as in Land et al. and Margara et al., and we tuned the same parameters (k_u, ntm, ν, µ and [Ca²⁺]_T50) within the same ranges, except for the Ca²⁺ sensitivity [Ca²⁺]_T50. In both Land et al. and Margara et al., [Ca²⁺]_T50 was optimized within [0.8, 0.9]. However, values in that range would have altered too much the CaT amplitude of BPS 2020. Land et al. already reported that such parameter “needs to be fit depending on the calcium transient used to drive the model,” as it is not consistent inter-species and also variable in their experiments on skinned human CMs. Therefore, we optimized [Ca²⁺]_T50 in the range [0.5, 0.6], which allowed us to keep the same CaT morphology and magnitude of the original BPS 2020. As we reported in Section 2, we followed a hybrid optimization approach based on genetic algorithm (Step 1, as in Land et al. and Margara et al.) to avoid being stuck in local minima, followed by a simplex (Step 2) on the sarcoplasmic Ca²⁺ fluxes again to keep the BPSLand Ca²⁺ handling the most similar to BPS2020. Already at this stage, the resulting model would have satisfactorily simulated the considered AP and Ta biomarkers. However, it lost the ability to simulate the inverse APD-[Ca²⁺]_o dependence for high [Ca²⁺]_o values. Such dependence was one of the key-novelties of BPS2020 (Bartolucci et al., 2020). In order to restore it (Figure 1), we added one step of manual tuning on the sarcoplasmic maximal fluxes J_rel,max and J_up,max, applying only minor changes fully consistent with the physiological formulation (Step 3). The final BPSLand model satisfactorily simulates AP and Ta biomarkers, together with the APD-[Ca²⁺]_o inverse dependence (Figure 1; Table 4).

4.2 Validation of the Model Against In Vitro Data and Comparison With Other In Silico Models

We first validated BPSLand against the same AP data (APD rate dependence and restitution in control condition and with current blocker) used to validate BPS2020 and the original ORd model. The rationale is we want BPSLand to work as well as BPS2020 in simulating electrophysiology data. As we already presented in detail those simulations in Bartolucci et al. (2020), here we report our results and the used protocols in the Supplementary Section S1 and Supplementary Figures S2, S3). These results confirm that adding the mechanical model have not altered the behaviour of the model electrophysiology. Nevertheless, it should be taken into account that having a single experimental dataset, including both electrical and mechanical measurements, would be the ideal setting to better calibrate and validate an electromechanical model (same in vitro preparations, clearer assessment of the mechanoelectric feedback, etc.). However, to our current knowledge, there is no such kind of data collection.

In terms of contractile properties, we compared BPSLand simulations to in vitro experiments performed on human CMs and cardiac preparations (Section 2.1). BPSLand successfully simulated the linear force-frequency dependence reported by Janssen and Periasamy (2007) (Figure 2). Such dependence was previously simulated by Lyon et al. (2020), although obtaining lower values of normalized force compared to BPSLand and to in vitro data in the range [1, 2.5] Hz (see Figure 2B in the original Lyon et al. paper). In terms of relaxation time, BPSLand optimally replicated the CaRT₅₀ data by Pieske et al., and very well the TaRT₅₀ by Pieske et al. (1995) and Janssen and Periasamy (2007). BPSLand TaRT₅₀ is lower at the slowest pacing rates (Figure 2). We ascribe this discrepancy to the TaRT₅₀ interval we used at 1 Hz during the model optimization: BPSLand and Janssen07 TaRT₅₀ are positioned at the opposite sides of such interval (purple line in Figure 2) while Pieske95 is out of bound. In terms of transmural heterogeneity (Figure 3), BPSLand simulations are in agreement with the in silico results of ToR-ORd+Land in terms of APD (M > ENDO > EPI) and TaPeak (M > EPI > ENDO) sequences, although the TaPeak values are greater in ToR-ORd+Land than in our model (M ∼ 60 kPa, EPI ∼ 40 kPa, ENDO ∼ 20 kPa). In fact, although we used the same TaPeak range as in Margara et al. (2021), i.e., [15, 25] kPa, BPSLand simulates a reference ENDO TaPeak equals to 15.6 kPa, which is more in line with the Haynes et al. (2014) in vitro values (Figure 3), especially for the small difference we observed in our ENDO vs. EPI TaPeak. As in Margara et al. (2021), we have tested heterogeneities in myofilament calcium sensitivities by acting on the baseline [Ca²⁺]_T50 value for the epicardial cell type (Figure 3E), showing that a small change of the [Ca²⁺]_T50 parameter replicates better the experiments. This result suggests that simulating transmural heterogeneity with electromechanical models may not only require re-calibration of the electrophysiological part but also of the mechanical part of the chosen model (Haynes et al., 2014).

Abnormalities in the ionic regulations of cardiomyocytes e.g., EADs and DADs can trigger the occurrence of a contractile irregularity in form of aftercontractions (Nguyen et al., 2017) the incidence of which has been reported in animal models of heart failure associated with arrhythmogenesis (Pogwizd et al., 2001). We observed that the endocardial BPSLand, as the original BPS2020, reacts to dofetilide and quinidine not producing EADs nor aftercontractions, but with an extreme prolongation of APD. This is not surprising, since we designed BPSLand carefully maintaining the electrophysiology of BPS2020. On the other hand, the M cell model reacted to both drugs with such abnormalities in electrophysiology and contractility. From the modelling point of view, it is not surprising: compared to the endocardial model, the M cell model has smaller G_Kr (thus smaller repolarization reserve), larger G_CaL (thus being prone to more significant reactivation of I_CaL during phase 3 of the AP) and J_rel,max (i.e., larger releases, also spontaneous, of Ca²⁺ from SR). From the in vitro point of view, the higher sensitivity of M cells to drugs affecting repolarizing ion currents compared to endocardial and epicardial was reported by Antzelevitch et al. (1999), with a panel of 13 drugs. Nonetheless, we did not observe a 1:1 correspondence between EADs and aftercontractions. We previously observed (see Figure 6C in the original BPS2020 paper) EADs triggered by two different mechanisms: I_CaL reactivation-driven and RyR spontaneous opening-driven EADs, as we also report here in Figure 5. Only in the case of a RyR spontaneous opening-driven EAD, we also have the corresponding aftercontraction, which is not present for I_CaL-driven EADs. Similarly to RyR spontaneous opening-driven EADs, also DADs are source of aftercontractions (Desantiago et al., 2008). BPSLand correctly simulated them using a protocol aimed to stress the model. Aftercontractions have been reported in vitro following the administration of diverse compounds or in presence of mutations in several cardiac preparations, e.g., cardiac tissues and trabeculae following dofetilide administration (Nguyen et al., 2017), in myocardial slices containing titin and collagen administered with isoproterenol (Watson et al., 2019), or in CPVT human induced pluripotent stem cell-derived CMs (Novak et al., 2012).

4.3 Limitations

The proposed computational model can be used to better understand the electromechanical interactions and the strong relationship between Ca²⁺ regulation and mechanics. Despite this, the experimental in vitro human data, taking into account both electrical and mechanical aspects, are still few, and urgently required to ensure better insight in electromechanical coupling and design more accurate models. The BPSLand model itself has some limitations. Preload and afterload conditions contribute to contractility response and should be considered in future model developments by including a mathematical description of dynamic changes in sarcomere length, since now only the isometric condition can be simulated. Other previously published mechanical models, e.g., Rice et al. (2008), Dupuis et al. (2016), and Dupuis et al. (2018), also include a mechanical description of sarcomere lengthening and shortening, thus expanding the range of possible simulations. Our choice to use Land model is based on the fact that it is validated against human experimental data. BPSLand model describes mechano-electric feedback only through the binding of calcium to troponin, but this phenomenon also includes other actors, for example stretch-activated ion channels (Peyronnet et al., 2016), which are modulated by membrane stretch and yield a current acting on the cardiomyocyte membrane potential. Future works should include into the model also these channels. Figures 2B,C show another limitation of BPSLand: while the model captures well the CaRT₅₀ in vitro data, it slightly underestimates the TaRT₅₀, as BPSLand simulates very similar CaRT₅₀ and TaRT₅₀ for each tested rate. One reason could be a slightly too fast relaxation dynamic in the contractile element. However, we replicated the same test with a second in silico model (Supplementary Figure S4) and an even more different behavior emerges. As cardiomyocyte electromechanical models are not so common yet as electrophysiology only models, it is clear that further iterations of optimization and validation shall be made in the future. Finally, we did not test the application of the model to multiscale simulations (2D or 3D) since it was beyond the scope of the work, although it will be interesting to check BPSLand behaviour also in this field of applicability.

5 Conclusion

In this paper, we presented our new electromechanical model of human adult ventricular cardiomyocyte, built and validated using several sets of human in vitro experiments. In addition to replicate correctly the results produced by its predecessor BPS2020, BPSLand adds an accurate simulation of active tension and contractility abnormalities, which can be triggered by drugs or specific pacing protocols. Therefore, BPSLand expands the domain of applicability of in silico model, which traditionally focus mainly on the simulation of the cardiac cell electrophysiology.

Data Availability Statement

The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

Author Contributions

CB and MP conceived and designed the study. CB and MP developed and validated the in silico model. MF, CB, SS, and MP analysed the in silico data, prepared the figures, and drafted the manuscript. All the authors interpreted the results and revised the manuscript.

Funding

CB was supported in part by the PRIN (Progetti di Ricerca di Rilevante Interesse Nazionale) 2017, Prot. N.2017AXL54F. MF was supported by the Graduate School of Faculty of Medicine and Health Technology, Tampere University. MP was supported by the Finnish Cultural Foundation, decision 210813, and by Academy of Finland Centre of Excellence in Body-on-Chip Research.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The authors thank the Tampere Center for Scientific Computing for generous computational resources.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2022.906146/full#supplementary-material

References

Ahola A., Pölönen R.-P., Aalto-Setälä K., Hyttinen J. (2018). Simultaneous Measurement of Contraction and Calcium Transients in Stem Cell Derived Cardiomyocytes. Ann. Biomed. Eng. 46, 148–158. doi:10.1007/s10439-017-1933-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Antzelevitch C., Shimizu W., Yan G.-X., Sicouri S., Weissenburger J., Nesterenko V. V. (1999). The M Cell: Its Contribution to the ECG and to Normal and Abnormal Electrical Function of the Heart. J. Cardiovasc. Electrophysiol. 10, 1124–1152. doi:10.1111/j.1540-8167.1999.tb00287.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Awinda P. O., Bishaw Y., Watanabe M., Guglin M. A., Campbell K. S., Tanner B. C. W. (2020). Effects of Mavacamten on Ca2+ Sensitivity of Contraction as Sarcomere Length Varied in Human Myocardium. Br. J. Pharmacol. 177, 5609–5621. doi:10.1111/bph.15271

PubMed Abstract | CrossRef Full Text | Google Scholar

Bartolucci C., Passini E., Hyttinen J., Paci M., Severi S. (2020). Simulation of the Effects of Extracellular Calcium Changes Leads to a Novel Computational Model of Human Ventricular Action Potential with a Revised Calcium Handling. Front. Physiol. 11, 1–20. doi:10.3389/fphys.2020.00314

PubMed Abstract | CrossRef Full Text | Google Scholar

Campbell S. G., Lionetti F. V., Campbell K. S., McCulloch A. D. (2010). Coupling of Adjacent Tropomyosins Enhances Cross-Bridge-Mediated Cooperative Activation in a Markov Model of the Cardiac Thin Filament. Biophys. J. 98, 2254–2264. doi:10.1016/j.bpj.2010.02.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Clancy C. E., Rudy Y. (2002). Na+ Channel Mutation that Causes Both Brugada and Long-QT Syndrome Phenotypes: A Simulation Study of Mechanism. Circulation 105, 1208–1213. doi:10.1161/hc1002.105183

PubMed Abstract | CrossRef Full Text | Google Scholar

D’Errico J. (2022). Fminsearchbnd, Fminsearchcon. MATLAB Cent. File Exch. Available at: https://www.mathworks.com/matlabcentral/fileexchange/8277-fminsearchbnd-fminsearchcon.

Google Scholar

Desantiago J., Ai X., Islam M., Acuna G., Ziolo M. T., Bers D. M., et al. (2008). Arrhythmogenic Effects of β2-adrenergic Stimulation in the Failing Heart Are Attributable to Enhanced Sarcoplasmic Reticulum Ca Load. Circ. Res. 102, 1389–1397. doi:10.1161/CIRCRESAHA.107.169011

PubMed Abstract | CrossRef Full Text | Google Scholar

Dupuis L. J., Lumens J., Arts T., Delhaas T. (2018). High Tension in Sarcomeres Hinders Myocardial Relaxation: A Computational Study. PLoS One 13, e0204642. doi:10.1371/journal.pone.0204642

PubMed Abstract | CrossRef Full Text | Google Scholar

Dupuis L. J., Lumens J., Arts T., Delhaas T. (2016). Mechano-chemical Interactions in Cardiac Sarcomere Contraction: A Computational Modeling Study. PLoS Comput. Biol. 12, 1–20. doi:10.1371/journal.pcbi.1005126

PubMed Abstract | CrossRef Full Text | Google Scholar

Dutta S., Chang K. C., Beattie K. A., Sheng J., Tran P. N., Wu W. W., et al. (2017). Optimization of an In Silico Cardiac Cell Model for Proarrhythmia Risk Assessment. Front. Physiol. 8, 1–15. doi:10.3389/fphys.2017.00616

PubMed Abstract | CrossRef Full Text | Google Scholar

Forouzandehmehr M., Ghoytasi I., Shamloo A., Ghosi S. (2022). Particles in Coronary Circulation: A Review on Modelling for Drug Carrier Design. Mat. Des. 216, 110511. doi:10.1016/j.matdes.2022.110511

CrossRef Full Text | Google Scholar

Fülöp G., Oláh A., Csipo T., Kovács Á., Pórszász R., Veress R., et al. (2021). Omecamtiv Mecarbil Evokes Diastolic Dysfunction and Leads to Periodic Electromechanical Alternans. Basic Res. Cardiol. 116, 1–16. doi:10.1007/s00395-021-00866-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo D., Liu Q., Liu T., Elliott G., Gingras M., Kowey P. R., et al. (2011). Electrophysiological Properties of HBI-3000: A New Antiarrhythmic Agent with Multiple-Channel Blocking Properties in Human Ventricular Myocytes. J. Cardiovasc. Pharmacol. 57, 79–85. doi:10.1097/FJC.0b013e3181ffe8b3

PubMed Abstract | CrossRef Full Text | Google Scholar

Haynes P., Nava K. E., Lawson B. A., Chung C. S., Mitov M. I., Campbell S. G., et al. (2014). Transmural Heterogeneity of Cellular Level Power Output Is Reduced in Human Heart Failure. J. Mol. Cell. Cardiol. 72, 1–8. doi:10.1016/j.yjmcc.2014.02.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Janssen P. M. L., Periasamy M. (2007). Determinants of Frequency-dependent Contraction and Relaxation of Mammalian Myocardium. J. Mol. Cell. Cardiol. 43, 523–531. doi:10.1016/j.yjmcc.2007.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Katz A. M. (1973). Effects of Ischemia on the Contractile Processes of Heart Muscle. Am. J. Cardiol. 32, 456–460. doi:10.1016/S0002-9149(73)80036-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Kernik D. C., Yang P.-C., Kurokawa J., Wu J. C., Clancy C. E. (2020). A Computational Model of Induced Pluripotent Stem-Cell Derived Cardiomyocytes for High Throughput Risk Stratification of KCNQ1 Genetic Variants. PLOS Comput. Biol. 16, e1008109. doi:10.1371/journal.pcbi.1008109

PubMed Abstract | CrossRef Full Text | Google Scholar

Koivumäki J. T., Naumenko N., Tuomainen T., Takalo J., Oksanen M., Puttonen K. A., et al. (2018). Structural Immaturity of Human iPSC-Derived Cardiomyocytes: In Silico Investigation of Effects on Function and Disease Modeling. Front. Physiol. 9, 1–17. doi:10.3389/fphys.2018.00080

PubMed Abstract | CrossRef Full Text | Google Scholar

Land S., Niederer S. A., Louch W. E., Sejersted O. M., Smith N. P. (2013). Integrating Multi-Scale Data to Create a Virtual Physiological Mouse Heart. Interface Focus 3. doi:10.1098/rsfs.2012.0076

PubMed Abstract | CrossRef Full Text | Google Scholar

Land S., Park-Holohan S. J., Smith N. P., dos Remedios C. G., Kentish J. C., Niederer S. A. (2017). A Model of Cardiac Contraction Based on Novel Measurements of Tension Development in Human Cardiomyocytes. J. Mol. Cell. Cardiol. 106, 68–83. doi:10.1016/j.yjmcc.2017.03.008

PubMed Abstract | CrossRef Full Text | Google Scholar

Leitch S. (1996). Effect of Raised Extracellular Calcium on Characteristics of the Guinea-pig Ventricular Action Potential. J. Mol. Cell. Cardiol. 28, 541–551. doi:10.1006/jmcc.1996.0050

PubMed Abstract | CrossRef Full Text | Google Scholar

Li P., Rudy Y. (2011). A Model of Canine Purkinje Cell Electrophysiology and Ca 2+ Cycling. Circ. Res. 109, 71–79. doi:10.1161/CIRCRESAHA.111.246512

PubMed Abstract | CrossRef Full Text | Google Scholar

Li Z., Mirams G. R., Yoshinaga T., Ridder B. J., Han X., Chen J. E., et al. (2020). General Principles for the Validation of Proarrhythmia Risk Prediction Models: An Extension of the CiPA In Silico Strategy. Clin. Pharmacol. Ther. 107, 102–111. doi:10.1002/cpt.1647

PubMed Abstract | CrossRef Full Text | Google Scholar

Lyon A., Dupuis L. J., Arts T., Crijns H. J. G. M., Prinzen F. W., Delhaas T., et al. (2020). Differentiating the Effects of β-adrenergic Stimulation and Stretch on Calcium and Force Dynamics Using a Novel Electromechanical Cardiomyocyte Model. Am. J. Physiol. - Hear. Circ. Physiol. 319, H519–H530. doi:10.1152/ajpheart.00275.2020

CrossRef Full Text | Google Scholar

Margara F., Wang Z. J., Levrero-Florencio F., Santiago A., Vázquez M., Bueno-Orovio A., et al. (2021). In-silico Human Electro-Mechanical Ventricular Modelling and Simulation for Drug-Induced Pro-arrhythmia and Inotropic Risk Assessment. Prog. Biophys. Mol. Biol. 159, 58–74. doi:10.1016/j.pbiomolbio.2020.06.007

PubMed Abstract | CrossRef Full Text | Google Scholar

Mulieri L. A., Hasenfuss G., Leavitt B., Allen P. D., Alpert N. R. (1992). Altered Myocardial Force-Frequency Relation in Human Heart Failure. Circulation 85, 1743–1750. doi:10.1161/01.CIR.85.5.1743

PubMed Abstract | CrossRef Full Text | Google Scholar

Musuamba F. T., Skottheim Rusten I., Lesage R., Russo G., Bursi R., Emili L., et al. (2021). Scientific and Regulatory Evaluation of Mechanistic In Silico Drug and Disease Models in Drug Development: Building Model Credibility. CPT Pharmacometrics Syst. Pharmacol. 10, 804–825. doi:10.1002/psp4.12669

PubMed Abstract | CrossRef Full Text | Google Scholar

Nagy N., Acsai K., Kormos A., Sebők Z., Farkas A. S., Jost N., et al. (2013). [Ca2+]i-induced Augmentation of the Inward Rectifier Potassium Current (Ik1) in Canine and Human Ventricular Myocardium. Pflügers Arch. - Eur. J. Physiol. 465, 1621–1635. doi:10.1007/s00424-013-1309-x

CrossRef Full Text | Google Scholar

Nguyen N., Nguyen W., Nguyenton B., Ratchada P., Page G., Miller P. E., et al. (2017). Adult Human Primary Cardiomyocyte-Based Model for the Simultaneous Prediction of Drug-Induced Inotropic and Pro-arrhythmia Risk. Front. Physiol. 8, 1–16. doi:10.3389/fphys.2017.01073

PubMed Abstract | CrossRef Full Text | Google Scholar

Niederer S. A., Lumens J., Trayanova N. A. (2019). Computational Models in Cardiology. Nat. Rev. Cardiol. 16, 100–111. doi:10.1038/s41569-018-0104-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Novak A., Barad L., Zeevi-Levin N., Shick R., Shtrichman R., Lorber A., et al. (2012). Cardiomyocytes Generated from CPVT D307H Patients Are Arrhythmogenic in Response to β-adrenergic Stimulation. J. Cell. Mol. Med. 16, 468–482. doi:10.1111/j.1582-4934.2011.01476.x

PubMed Abstract | CrossRef Full Text | Google Scholar

O’Hara T., Virág L., Varró A., Rudy Y. (2011). Simulation of the Undiseased Human Cardiac Ventricular Action Potential: Model Formulation and Experimental Validation. PLoS Comput. Biol. 7, e1002061. doi:10.1371/journal.pcbi.1002061

PubMed Abstract | CrossRef Full Text | Google Scholar

Paci M., Casini S., Bellin M., Hyttinen J., Severi S. (2018a). Large-Scale Simulation of the Phenotypical Variability Induced by Loss-Of-Function Long QT Mutations in Human Induced Pluripotent Stem Cell Cardiomyocytes. Int. J. Mol. Sci. 19, 3583. doi:10.3390/ijms19113583

PubMed Abstract | CrossRef Full Text | Google Scholar

Paci M., Koivumäki J. T., Lu H. R., Gallacher D. J., Passini E., Rodriguez B. (2021). Comparison of the Simulated Response of Three In Silico Human Stem Cell-Derived Cardiomyocytes Models and In Vitro Data under 15 Drug Actions. Front. Pharmacol. 12, 1–16. doi:10.3389/fphar.2021.604713

CrossRef Full Text | Google Scholar

Paci M., Passini E., Klimas A., Severi S., Hyttinen J., Rodriguez B., et al. (2018b). “In Silico Populations Optimized on Optogenetic Recordings Predict Drug Effects in Human Induced Pluripotent Stem Cell-Derived Cardiomyocytes,” in 2018 Computing in Cardiology Conference (CinC),Computing in Cardiology Conference (CinC). doi:10.22489/CinC.2018.086

CrossRef Full Text | Google Scholar

Paci M., Passini E., Severi S., Hyttinen J., Rodriguez B. (2017). Phenotypic Variability in LQT3 Human Induced Pluripotent Stem Cell-Derived Cardiomyocytes and Their Response to Antiarrhythmic Pharmacologic Therapy: An In Silico Approach. Hear. Rhythm 14, 1704–1712. doi:10.1016/j.hrthm.2017.07.026

CrossRef Full Text | Google Scholar

Papadaki M., Kampaengsri T., Barrick S. K., Campbell S. G., von Lewinski D., Rainer P. P., et al. (2022). Myofilament Glycation in Diabetes Reduces Contractility by Inhibiting Tropomyosin Movement, Is Rescued by cMyBPC Domains. J. Mol. Cell. Cardiol. 162, 1–9. doi:10.1016/j.yjmcc.2021.08.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Passini E., Britton O. J., Lu H. R., Rohrbacher J., Hermans A. N., Gallacher D. J., et al. (2017). Human In Silico Drug Trials Demonstrate Higher Accuracy Than Animal Models in Predicting Clinical Pro-arrhythmic Cardiotoxicity. Front. Physiol. 8, 1–15. doi:10.3389/fphys.2017.00668

PubMed Abstract | CrossRef Full Text | Google Scholar

Peyronnet R., Nerbonne J. M., Kohl P. (2016). Cardiac Mechano-Gated Ion Channels and Arrhythmias. Circ. Res. 118, 311–329. doi:10.1161/CIRCRESAHA.115.305043

PubMed Abstract | CrossRef Full Text | Google Scholar

Pieske B., Kretschmann B., Meyer M., Holubarsch C., Weirich J., Posival H., et al. (1995). Alterations in Intracellular Calcium Handling Associated with the Inverse Force-Frequency Relation in Human Dilated Cardiomyopathy. Circulation 92, 1169–1178. doi:10.1161/01.CIR.92.5.1169

PubMed Abstract | CrossRef Full Text | Google Scholar

Pieske B., Sütterlin M., Schmidt-Schweda S., Minami K., Meyer M., Olschewski M., et al. (1996). Diminished Post-rest Potentiation of Contractile Force in Human Dilated Cardiomyopathy: Functional Evidence for Alterations in Intracellular Ca2+ Handling. J. Clin. Invest. 98, 764–776. doi:10.1172/JCI118849

PubMed Abstract | CrossRef Full Text | Google Scholar

Pogwizd S. M., Schlotthauer K., Li L., Yuan W., Bers D. M. (2001). Arrhythmogenesis and Contractile Dysfunction in Heart Failure. Circ. Res. 88, 1159–1167. doi:10.1161/hh1101.091193

PubMed Abstract | CrossRef Full Text | Google Scholar

Rahman M. A., Ušaj M., Rassier D. E., Månsson A. (2018). Blebbistatin Effects Expose Hidden Secrets in the Force-Generating Cycle of Actin and Myosin. Biophys. J. 115, 386–397. doi:10.1016/j.bpj.2018.05.037

PubMed Abstract | CrossRef Full Text | Google Scholar

Rice J. J., Wang F., Bers D. M., De Tombe P. P. (2008). Approximate Model of Cooperative Activation and Crossbridge Cycling in Cardiac Muscle Using Ordinary Differential Equations. Biophys. J. 95, 2368–2390. doi:10.1529/biophysj.107.119487

PubMed Abstract | CrossRef Full Text | Google Scholar

Rossman E. I., Petre R. E., Chaudhary K. W., Piacentino V., Janssen P. M. L., Gaughan J. P., et al. (2004). Abnormal Frequency-dependent Responses Represent the Pathophysiologic Signature of Contractile Failure in Human Myocardium. J. Mol. Cell. Cardiol. 36, 33–42. doi:10.1016/j.yjmcc.2003.09.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Santini L., Palandri C., Nediani C., Cerbai E., Coppini R. (2020). Modelling Genetic Diseases for Drug Development: Hypertrophic Cardiomyopathy. Pharmacol. Res. 160, 105176. doi:10.1016/j.phrs.2020.105176

PubMed Abstract | CrossRef Full Text | Google Scholar

Sarkar S. S., Trivedi D. V., Morck M. M., Adhikari A. S., Pasha S. N., Ruppel K. M., et al. (2020). The Hypertrophic Cardiomyopathy Mutations R403Q and R663H Increase the Number of Myosin Heads Available to Interact with Actin. Sci. Adv. 6, 1–10. doi:10.1126/sciadv.aax0069

CrossRef Full Text | Google Scholar

Severi S., Corsi C., Cerbai E. (2009). From In Vivo Plasma Composition to In Vitro Cardiac Electrophysiology and In Silico Virtual Heart: the Extracellular Calcium Enigma. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 367, 2203–2223. doi:10.1098/rsta.2009.0032

CrossRef Full Text | Google Scholar

Sheikh F., Ouyang K., Campbell S. G., Lyon R. C., Chuang J., Fitzsimons D., et al. (2012). Mouse and Computational Models Link Mlc2v Dephosphorylation to Altered Myosin Kinetics in Early Cardiac Disease. J. Clin. Invest. 122, 1209–1221. doi:10.1172/JCI61134

PubMed Abstract | CrossRef Full Text | Google Scholar

Toepfer C. N., Garfinkel A. C., Venturini G., Wakimoto H., Repetti G., Alamo L., et al. (2020). Myosin Sequestration Regulates Sarcomere Function, Cardiomyocyte Energetics, and Metabolism, Informing the Pathogenesis of Hypertrophic Cardiomyopathy. Circulation 141, 828–842. doi:10.1161/CIRCULATIONAHA.119.042339

PubMed Abstract | CrossRef Full Text | Google Scholar

Tomek J., Bueno-Orovio A., Passini E., Zhou X., Minchole A., Britton O., et al. (2019). Development, Calibration, and Validation of a Novel Human Ventricular Myocyte Model in Health, Disease, and Drug Block. Elife 8, 1–48. doi:10.7554/eLife.48890

CrossRef Full Text | Google Scholar

Watson S. A., Duff J., Bardi I., Zabielska M., Atanur S. S., Jabbour R. J., et al. (2019). Biomimetic Electromechanical Stimulation to Maintain Adult Myocardial Slices In Vitro. Nat. Commun. 10. doi:10.1038/s41467-019-10175-3

PubMed Abstract | CrossRef Full Text | Google Scholar