Altered contractility in mutation-specific hypertrophic cardiomyopathy: A mechano-energetic in silico study with pharmacological insights

Abstract

Introduction: Mavacamten (MAVA), Blebbistatin (BLEB), and Omecamtiv mecarbil (OM) are promising drugs directly targeting sarcomere dynamics, with demonstrated efficacy against hypertrophic cardiomyopathy (HCM) in (pre)clinical trials. However, the molecular mechanism affecting cardiac contractility regulation, and the diseased cell mechano-energetics are not fully understood yet.

Methods: We present a new metabolite-sensitive computational model of human induced pluripotent stem cell-derived cardiomyocytes (hiPSC-CMs) electromechanics to investigate the pathology of R403Q HCM mutation and the effect of MAVA, BLEB, and OM on the cell mechano-energetics.

Results: We offer a mechano-energetic HCM calibration of the model, capturing the prolonged contractile relaxation due to R403Q mutation (∼33%), without assuming any further modifications such as an additional Ca²⁺ flux to the thin filaments. The HCM model variant correctly predicts the negligible alteration in ATPase activity in R403Q HCM condition compared to normal hiPSC-CMs. The simulated inotropic effects of MAVA, OM, and BLEB, along with the ATPase activities in the control and HCM model variant agree with in vitro results from different labs. The proposed model recapitulates the tension-Ca²⁺ relationship and action potential duration change due to 1 µM OM and 5 µM BLEB, consistently with in vitro data. Finally, our model replicates the experimental dose-dependent effect of OM and BLEB on the normalized isometric tension.

Conclusion: This work is a step toward deep-phenotyping the mutation-specific HCM pathophysiology, manifesting as altered interfilament kinetics. Accordingly, the modeling efforts lend original insights into the MAVA, BLEB, and OM contributions to a new interfilament balance resulting in a cardioprotective effect.

Introduction

With a prevalence ranging from of 1/500 to 1/200 (Semsarian et al., 2015; Malhotra and Sharma, 2017), HCM represents the most prevalent genetic cardiac disorder mainly associated with pathogenic variants in sarcomere protein genes (Santini et al., 2020). Pathologies such as myocardium hypercontractility (Sarkar et al., 2020), impaired relaxation (Toepfer et al., 2020), elevated cardiac energy consumption, diastolic dysfunction, arrhythmogenesis, and heart failure (Sarkar et al., 2020) manifest due to such variants. The driver of the cyclic interactions between thin and thick filaments is the ATP hydrolysis by myosin-the enzymatic motor of sarcomere (Spudich, 2014). In HCM, myosin binding protein C and adult cardiac myosin isoforms (primarily encoded by MYBPS3 and MYH7 genes, respectively) host most of these pathogenic variants in sarcomere (Ho et al., 2018; Toepfer et al., 2019; Schmid and Toepfer, 2021).

Mavacamten (MAVA), Blebbistatin (BLEB), and Omecamtiv mecarbil (OM) are compounds directly modulating myofilament dynamics with promising effectiveness in treatment of sarcomeric cardiomyopathies. MAVA is an allosteric inhibitor of cardiac myosin ATPase with a negative inotropic effect and demonstrating efficacy in R403Q HCM clinical trials (Green et al., 2016; Santini et al., 2020). BLEB, a well characterized ATPase inhibitor, alters the Ca²⁺ sensitivity of the myofilament and has been widely used in trials (Kampourakis et al., 2018; Rahman et al., 2018; Wang et al., 2018; Gyimesi et al., 2021). OM is a recently developed myosin ATPase activator with a positive inotropic effect enhancing cardiac contractility (Tsukamoto, 2019).

The mechanisms of action of these drugs and their effects on cardiomyocyte electro-mechano-energetics are still to be fully known and under active research (Tsukamoto, 2019; Santini et al., 2020). Accordingly, computational studies on the effects of these drugs are mostly lacking. Focusing on R403Q HCM mutation, Margara et al. (2021) investigated the efficacy of MAVA, simulating the tension-Ca²⁺ relationship and active tension curves. They hypothesized that the impaired tension relaxation phase is caused by feedback from crossbridge (XB) cycling to the thin filament. Although this assumption results in consistent simulated impaired tension relaxation, the proposed MAVA mechanism of action coupled with a lack of a metabolite-sensitive mechanism in the contractile element (CE) urged us to investigate further the cause of this impairment. On the other hand, in silico models have been employed to probe the effect of BLEB beside animal muscle fibre experiments (Rahman et al., 2018), highlighting the role of BLEB in shifting the rate-limiting momentum from weakly to strongly bound states in XB cycling (Rahman et al., 2018). Finally, as an ion channel study, the effect of OM has been simulated using an in silico model of human ventricular action potential (AP) focusing on pro-arrhythmic assessments (Qu et al., 2021). Of note, Qu et al. report an IC₅₀ of 125.5 µM highlighting the significance of OM influence on the CE in 1–10 µM of OM compared with channel blocking formalism (Qu et al., 2021). In summary, to the best of our knowledge, no computational study has reported experimentally validated drug-induced Ca²⁺ sensitivity, ATPase dynamics, and dose-dependent effect of MAVA, BLEB, and OM on the tension-Ca²⁺ relationship, regarding the HCM pathophysiology. This coupled with the facts that these compounds are essentially myosin ATPase activators/inhibitors necessitate studying the drug effects using advanced metabolite-sensitive models. In addition, HCM mutations misregulate sarcomere function and cardiac energy consumption (van der Velden et al., 2018). This further highlights the need for a model that enables capturing the electro-mechano-energetics in pathophysiological investigations.

Here, we investigate the pathophysiology of HCM R403Q myosin mutation using a computational metabolite-sensitive model of hiPSC-CMs electromechanics developed based on ordinary differential equations (ODEs). This model, named hiMCE, is an update of our previous model of hiPSC-CM electromechanics (Forouzandehmehr et al., 2021) capturing ATPase activity and accounting for metabolite-sensitive kinetics (Figure 1) in the XB cycling and extending the capacity of models of the molecular mechanism of contraction and the drug effect predictions. The drugs studied here are modulators of myofilament dynamics, we reparametrized the CE of the model using available experimental data and presented novel mechanistic methods to simulate the effect of MAVA, OM, and BLEB on the Ca²⁺ sensitivity, contractility, and energetics of the hiPSC-CMs.

FIGURE 1

Methods

Extension to a metabolite-sensitive contractile element

We previously integrated a reparametrized mathematical model of the CE by Rice et al. (2008) with a new passive force handling into the hiPSC-CM model of electrophysiology by Paci et al. (2020) and studied the inotropic effect of different compounds (Forouzandehmehr et al., 2021). Based on (Rice et al. (2008) CE, Tran et al. (2017) introduced a mathematical model of a CE that incorporated metabolic-sensitivity. This was achieved by extending the model by Rice et al. (2008), with new parameters that account for the competitive binding of metabolic protons (H⁺) to the binding sites of Ca²⁺ on troponin C, and incorporates the binding kinetics of MgADP in the XB cycling (Figure 1). The extended mechanistic description divides the strongly bound state post isomerized rotation (XB_pstR) into two substates in rapid equilibrium, AM1 and AM2, to capture MgADP binding kinetics in the XBs (Figure 1).

The thermodynamically constrained model of XB kinetics and Ca²⁺ activation is divided into four states including a non-permissive (N_XB), a permissive (P_XB), a pre power stroke state (XB_preR), and a post power stroke state (XB_pstR) (Figure 1) (Tran et al., 2017). XB_preR and XB_pstR states both denote the state of strongly bound myosin heads to the actins (Rice et al., 2008; Tran et al., 2017). In diastole, XBs are in N_XB state and when activated by Ca²⁺ XBs enter P_XB state where they can participate in the binding and unbinding of myosin heads and processes of tension generation (Figure 1) (Tran et al., 2017). The active tension produced as the result of the process is equal to the product of myosin head strain and strongly bound states fractional occupancy (Tran et al., 2017).

In the present work, we have taken this metabolite-sensitive CE model proposed by Tran et al. (2017), and modified it to incorporate active contraction mechanisms integrated with the Paci2020 model of hiPSC-CMs electrophysiology (Paci et al., 2020). This extended CE model was manually tuned using the information from a previous sensitivity analysis of the contractile machinery (Forouzandehmehr et al., 2021) (Table 1). We calibrated the model to capture the AP, CaT, and contractile experimental biomarkers measured in hiPSC-CMs in control conditions, listed in Table 6.

This calibrated control/healthy model variant was then modified to describe an HCM mutation and the effect of pharmaceutical mechanical modulators, and in sections 2.2 and 2.3, we explain the tuning of the relevant CE parameters for each scenario. To enlighten the role of these parameters, we describe the XB kinetics (Tran et al., 2010) as follows in relation to Figure 1: and are transition rates between and and responsible for XB cycling activation (Rice et al., 2008). takes into account the Pi-dependent transition rate in the XB and has been defined in (Tran et al., 2010). is equal to the attachment rate to (Tran et al., 2010) which is given as:

value is set to 500 s⁻¹ and is the temperature-dependence set to 6.25 (Rice et al., 2008), in this work. (species-dependent XB cycling rate scaler (Rice et al., 2008)) value is given in Table 1. denotes a thermodynamically constrained transition rate from to state accounting for MgATP release when the myosin heads transit from weakly to strongly attached XB states detailed in (Tran et al., 2010). is the MgADP and MgATP-sensitive transition rate from to (Tran et al., 2010) and is calculated as:Where denotes the reference (physiological) concentration of MgADP of 36 µM and is the MgADP dissociation constant detailed in (Tran et al., 2010). denotes a first order rate constant tuned by Tran et al. (2010) to maintain the validity of original Rice CE model (Rice et al., 2008) under physiological metabolic conditions regarding MgADP kinetics detailed in. is equal to the transition from to which is proton and MgADP-sensitive (Tran et al., 2010):Again, represents the adjustment of the rate constant affecting the transition between and states to include physiological proton and MgADP dependent effects further detailed in (Tran et al., 2010). Following the method proposed in (Margara et al., 2021), we defined F1 and F2 (Eq. (2)) as modulators of the transition between permissive binding state on actin (P_XB) and strongly bound XBs before isomerized rotation (XB_preR). In our model, F1 and F2 represent the (DRX:SRX)/(DRX:SRX)_control ratio (=1 in control mode) and indirectly affect the ATPase dynamics (Figure 1). is equal to the forward transition rate between XB_preR and XB_pstR, which is defined as (Rice et al., 2008):

The value for is given in Table 1. represents the temperature dependence, here set to 6.25. incorporates strain dependence into the forward transition rate, . is set to five and specifies the extent to which the isomerization rate is influenced by the mean strain of XB_preR (). is the mean strain (distortion) of XB_pstR state when the net motion between actin and myosin filaments is absent. Here, is set to 0.007 µm.

Furthermore, to and to also affect the time rate of change in the mean strains, Eqs (10), (11), and the steady-state population of strongly bound XB states Eqs (12)–(14) (Tran et al., 2010):Where denotes the sarcomere length velocity and represents an empirical scaler equal to 2 (Rice et al., 2008). For the detailed explanation of the metabolite-sensitive XB cycling we refer the readers to (Tran et al., 2010).

The contractile element calibration for HCM model variant

Our baseline CE inherits the main effects of contractile metabolic products such as MgATP, MgADP, inorganic Phosphate (Pi), and H⁺ on the tension development mechanism from the original Tran et al. model (Tran et al., 2010; Tran et al., 2017). We used this baseline model to develop an HCM mutant variant (R403Q) model, with altered myofilament kinetics. This model variant was created by modifying specific metabolic parameters to achieve a simulated state consistent with experimental reports of R403Q HCM (Table 2).

To obtain the HCM model variant, we changed F1 and F2 values following (Margara et al., 2021) and in line with the sensitivity test given in Supplementary Figure S4. We also increased the value of MgADP concentration and the reference value of Pi. Finally, we changed the ap2 coefficient regarding the model sensitivity given in Supplementary Figure S2.

The contractile element calibration for drug-induced effects

Previously, (Margara et al. (2021) assumed in their computational study that MAVA mainly influences the transitions between XB_preR and P_XB states following the DRX:SRX disturbing theory reported in (Toepfer et al., 2020). This was implemented by introducing F1 and F2 (with default values of 1) coefficients (Eq. (2)) representing DRX:SRX ratios to the time-dependent description of XB_preR state.

In addition, we simulated the effect of MAVA not only by altering the values of F1 and F2 but also modifying the parameters listed in Table 1 to obtain a comprehensive and accurate simulation of the effect of 0.5 µM MAVA on our HCM R403Q model variant as given in Table 3. These further modifications, done as a manual parameter tuning, are based on the reported effect of MAVA on Ca²⁺ activation and binding process (Awinda et al., 2020), Pi (Alsulami and Marston, 2020), and ATPase activity (Green et al., 2016). Specifically, changes in K_on and n_perm values have been made regarding the sensitivity analyses given in (Forouzandehmehr et al., 2021). The F1 and F2 values were changed according to model sensitivity behavior given in Supplementary Figure S4. Further, and coefficients and A-E values were obtained by trial and error.

To indicate, A modulates P_XB to XB_preR transition, B influences Pi-dependent XB_preR to P_XB transition, C takes effect on XB_preR to XB_pstR transition, D affects proton-dependent XB_pstR to XB_preR transition, and E controls MgATP-dependent transition from XB_pstR to P_XB states. Correspondingly, the proposed modulations (A to E in Table 3) to Eqs (15)–(19), are in line with a disturbed interfilament signaling that affects the force-producing states of XB suggested in the etiology of R403Q (Nag et al., 2015). Our motivation for the A-E coefficient modifications was the role of xbmodsp parameter in contraction relaxation time observed in the sensitivity study before (Forouzandehmehr et al., 2021). As changing xbmodsp solely could not lead to an accurate simulation of impaired relaxation restored by 0.5 µM MAVA, we used A-E values to optimize the distribution of xbmodsp effect on the XB cycling.

Values of thin filament regulation and XB cycling parameters , , and temperature dependences , , were directly taken from (Rice et al., 2008). Also, is a strain-dependent rate modifier defined in (Rice et al., 2008).

The values of parameters changed in the model to simulate the effects of 5 µM BLEB and 1 µM OM are given in Table 4. F1 and F2 values were obtained with attention to sensitivity plots given in Supplementary Figure S4. Similarly, and coefficients were found regarding the model behavior shown in Supplementary Figures S2, S3, and S5. Lastly, values of parameters listed in Table 4 rows 7 to 15 were obtained manually according to the sensitivity reports given in (Forouzandehmehr et al., 2021).

Moreover, to capture the dose-dependent effect of BLEB and OM on the normalized tension, based on model sensitivity tests (Supplementary Figures S2-S5) and previous sensitivity analyses (Forouzandehmehr et al., 2021), we identified the main variable of the CE governing the maximum developed tension, , and identified coefficient values accordingly (Supplementary Tables S1 and S2). Supplementary Figures S2-S5 show the sensitivity of tension-Ca²⁺ relationships, active tensions, and ATPase rates to , , and R (=F1 = F2) coefficients. Our HCM and drug-induced calibrations are unique sets obtained from an informed manual tuning, with e.g., 10%, 20%, 50% increments, started with a first guess within assumed boundaries (e.g., 1/100 to 5.4 times of baseline values). The examples of the results showing the increments and the model response are given in Supplementary Figures S2-S5. Overall, these sensitivity analyses combined with our previous sensitivity investigations on the CE of our previous model (Forouzandehmehr et al., 2021) were the basis of the informed manual parameter tunings done to obtain HCM- and drug-induced calibrations in this work.

The experimental data for calibrations and validations

Table 5 gives the experimental data using which the results of this work have been calibrated and validated.

Results

The metabolite-sensitive model of hiPSC-CMs

All the equations were solved using ode15 s integrator of MATLAB with a maximum step size of 10^–3 and an initial step size of 2 × 10^–5. To reach the steady state, the results were obtained after 800 beats paced at 1 Hz unless otherwise mentioned.

First, we show that our computational model can correctly simulate the main AP, Ca²⁺ transients (CaTs), and active tension (AT) biomarkers as had been simulated by our previous electromechanical hiPSC-CM-CE model (Forouzandehmehr et al., 2021). As Table 6 shows, our new metabolite-sensitive hiMCE model is able to simulate the main biomarkers within the experimental ranges in the validation datasets. The increased thermodynamic detail of the CE did not significantly alter the biomarker values compared to our previous reparameterization (Forouzandehmehr et al., 2021) and the original Paci et al. hiPSC-CM model (Paci et al., 2020). Also, Figure 2 shows the contractility characteristics simulated using the hiMCE model are consistent with the previously validated results, while also illustrating selected fundamental outputs (Ruan et al., 2016; Pioner et al., 2020).

FIGURE 2

Pioner et al., 2020; Forouzandehmehr et al., 2021).

Hypercontractility in R403Q HCM and mavacamten

In order to simulate the abnormal prolonged relaxation in the developed active tension due to R403Q HCM mutation, Margara et al. (2021) hypothesized a feedback from XB cycling to the thin filament activation. To investigate further, using the parameter values in Table 2 and consistent with the metabolic data detailed in section 2.2, we simulated the active tension and ATPase rate in R403Q HCM model variant. The CaT morphology remains unchanged in the HCM R403Q mode (Figure 3B), consistently with experimental data reported for hiPSC-CMs in (Sewanan et al., 2019). Interestingly, the results in Figure 3C suggest that including energetics in the CE reacts to the pathological changes due to HCM and can correctly predict the prolonged relaxation in the developed active tension (∼33%), consistently with in vitro hiPSC-CMs data (Toepfer et al., 2020). Moreover, the increased fractional cell shortening (∼40%) due to the R403Q mutation is consistent with experimental measurements in (Toepfer et al., 2020). Of note, the model also correctly predicts the negligible change in the ATPase activity (Figure 3D), consistently with the experimental data (Table 5) (Nag et al., 2015; Sarkar et al., 2020).

FIGURE 3

Gollapudi et al., 2021).

To simulate the electro-mechano-energetic effect of 0.5 µM MAVA, we used the model calibration values listed in Table 3. We have assumed that MAVA would shift the elevated metabolites in the HCM model variant, Pi and MgADP, towards their baseline values. Our model could accurately predict the unaffected CaTs due to MAVA as reported experimentally earlier (Green et al., 2016). Further, the order of reduction in the simulated ATPase rate (19.3%) due to 0.5 µM MAVA, Figures 3D, 3is within the reduction range, 17.9–28.5%, reported in previous experimental ATPase activity measurements (Gollapudi et al., 2021). Also, the model consistently predicts the reduction in the relaxation phase in the ATPase rate (Figure 3D) (Kawas et al., 2017). Notably, our simulations quantitatively capture the reduction in the fractional cell shortening and prolonged tension relaxation due to R403Q mutated hiPSC-CMs after 0.5 µM MAVA (14.6% and 20.9%, respectively), consistently with recent experimental measurements (Toepfer et al., 2020). Finally as shown in Supplementary Figure S1, the CE model accurately predicts the unchanged pCa₅₀ in the tension-Ca²⁺ relationship consistent with the experimental data for 0.5 µM MAVA (Green et al., 2016).

Simulated effects of omecamtiv mecarbil and blebbistatin

Using the CE parameter values listed in Table 4, we simulated the effect of 5 µM BLEB and 1 µM OM. As Table 7 shows, the drug-induced calibration of the hiMCE model results in accurate predictions of Ca²⁺ sensitive effects of 5 µM BLEB and 1 µM OM consistent with experimental data (Kampourakis et al., 2018) as also Figure 4A qualitatively confirms. Additionally, the selected values for the coefficients of the tension governing variables in the CE, (Figures 4B,C and Supplementary Tables S1, S2), leads to the correct dose-dependent prediction for BLEB and OM (Figures 4D,E) and the expected inverse Hill curves reported experimentally (Kampourakis et al., 2018).

FIGURE 4

Kampourakis et al., 2018). T₀: isometric force in the absence of drugs.

Moreover, our model predicts an insignificant reduction in AP duration (3.4%) due to 1 µM OM (Figure 5A), evaluated by calculating APD₉₀ values (Figure 5A). This translates to 17 ms reduction in APD₉₀ (502–485 ms) which is consistent with the order of APD₉₀ reduction due to 1 µM OM, 12.2 ms, reported for canine cardiomyocytes at 1 Hz pacing (Szentandrassy et al., 2016). Also, the simulated increase in the basal ATPase rate due to 1 µM OM (Figure 5H) is qualitatively consistent with the experimental data reported before (Bakkehaug et al., 2015; Utter et al., 2015).

FIGURE 5

Our model predicts a 23% increase in the amplitude of the Ca²⁺ flux towards the myofilament, JCB (Figure 5B). The subsequent accumulation of intracellular Ca²⁺ is seen as a 4.5% increase in CaT peak (Figure 5E). Interestingly, the steady-state alterations in sarcolemmal Ca2+ transport are very subtle: virtually unchanged I_CaL (Figure 5C) and only very slightly increased I_pCa (Figure 5F). Whereas there is a more substantial 17% increase in the amplitude of the I_NCX (Figure 5D). The enhanced reverse mode of I_NCX causes accumulation of intracellular Na⁺ (7.04 vs. 7.43 mM) that promotes a stronger repolarizing I_NaK (Figure 5G). This appears to be the mechanism that causes the subtle yet visible 3.4% decrease in the AP duration (Figure 5A), consistently with the reported experiments suggesting OM as a safe compound on cardiac electrophysiology in clinically tolerated doses (Szentandrassy et al., 2016). These predictions can be insightful regarding the consequences of the disturbed interfilament signaling that OM elicits in the XBs.

Discussion

HCM and energetics of contraction

Cardiomyocytes, with no self-renewal capacity, must provide two billion beats during an average lifetime for which the cardiac muscle requires a significant amount of energy, 6 kg of ATP per day (Ingwall, 2002). This energy consumption is predominantly due to the function of sarcomeres in contraction. Therefore, the pathological conditions directly caused by sarcomeric mutations, such as R403Q HCM, necessitate studying contractile function of cardiomyocytes regarding cardiac metabolism. Our analysis demonstrates that the incorporated scheme of the metabolite-sensitive CE is able to capture the impaired (prolonged) tension relaxation, ∼33%, due to the R403Q mutation. Interestingly, with the energetics included, the additional feedback from XB cycling to the thin filaments, proposed previously by Margara et al. (2021), was not necessary to replicate the altered relaxation. This further highlights the importance of considering (patho)physiologically constrained metabolite-sensitive computational models in the investigation of sarcomeric cardiomyopathies.

HCM and drug-induced model calibrations

Contractile energetics become highly important when studying promising drugs reported in HCM clinical trials such as MAVA, BLEB, and OM. As our results show, the quantitatively valid simulation of the effect of MAVA, BLEB, and OM, in single dose or dose-dependently, could not be done without the calibration of parameters in the CE that directly or indirectly affect the energetics (Tables 3, 4, S1, and S2). Markedly, one of the important insights of this study stems from the parameters involved in the calibration of the model for the simulation of 0.5 µM MAVA. We took the MAVA modeling one step further by calibrating the CE altering parameters affecting Pi-dependent transition between permissive binding state on actin (P_XB) and the strongly bound XBs before isomerized rotation (XB_preR) state. Also, we modulated Ca²⁺ binding and sensitivity of the CE, and MgATP-dependent transitions between P_XB and XBs in strongly bound post isomerized rotation state (XB_postR) in accord with experimental metabolic reports (Green et al., 2016; Alsulami and Marston, 2020; Awinda et al., 2020). Towards decoding the precise drug mechanism of action, the modulations proposed here to explain the effect of MAVA (Table 3) implies that, alongside altering the disturbed DRX:SRX ratio, MAVA might also induce a new interfilament equilibrium, modulating tension-producing and energetic terms that explicitly affect the reverse transition at play between XB_preR and XB_postR affecting the strain-dependent isomerization of myosin heads. This possible pharmacological insight emerging from our model is interesting as MAVA mechanism of action inherently shifts the R403Q impaired metabolism towards normal regulation and this involves the impaired proton-dependent transition in R403Q mode (Toepfer et al., 2020).

The drug- and HCM-related calibrations presented in this work are in line with the proposed OM and BLEB structure-function relationships detailed in (Kampourakis et al., 2018). To enumerate, Kampourakis et al. (2018) have implied that OM-bound myosin heads relocate the tropomyosin to its on state when binding to actin in the absence of Ca²⁺ bound to troponin. Importantly, the XB activation, which is due to the effect of OM, has been significantly attributed to the stabilizing of the ON state of thick filaments. Further, these stabilized ON positions in thick filaments have been considered to promote an ON thin filament state through preventing tropomyosin returning to their off positions. This has been translated in our model by modifying k_off constants which are the rate constants affecting Ca²⁺ unbinding from low and high affinity sites on troponin. These transition rates affect regulatory sites on cardiac troponin leading to activation of XB cycle (Tran et al., 2010).

On the other hand, at intermediate Ca²⁺ in the physiological range (pCa nine to 4.3), OM activates actins along with the thick filaments. In our model, this has been translated by modifying K_on, K_np and K_pn rate constants for Ca²⁺ binding to troponin, forward and backward transition rates between P_xb and N_xb states, respectively (Table 4). Further, the constant rates directly affecting Ca²⁺ bound troponin thin filaments regulation induced by Ca²⁺ bound troponin (K_on, n_perm, k_offL, k_offH, k_pn, k_np, and k_offmode) have only been calibrated towards activation for OM and they are missing in BLEB calibration as BLEB does not switch both filaments to their ON states (Kampourakis et al., 2018). Further, modification of coefficient of Tropreg variable, the fraction of actins with Ca²⁺ bound, detailed in (Rice et al., 2008), in both OM and BLEB calibrations is consistent with the reported effect of these drugs as both OM and BLEB generally greatly decrease the Ca²⁺ activation co-operativity (Kampourakis et al., 2018).

In the HCM variant model, the significant increase in Pi concentration is consistent with the experimental metabolic reports of HCM R403Q and R92Q mutations in animal mouse models (Spindler et al., 1998; Abraham et al., 2013). Moreover, the impaired coronary perfusion due to HCM has been related to an abnormal energy reproduction that contributes to elevation of ADP and Pi (Sequeira et al., 2019), which is also consistent with conservation of phosphate and creatin reaction (Tran et al., 2010). Congruently, the HCM mutation-induced alterations in myofilament kinetics lead to increase in ADP-mediated products (van der Velden et al., 2018). Therefore, we increased the MgADP concentration within physiological ranges (Tran et al., 2010). Notably, the Ca²⁺ dependence in the activation of myofilaments for HCM and dilated cardiomyopathy has been shown to be altered in a similar fashion due to OM and BLEB influence (Spudich, 2014), implying that the etiology of HCM includes a disturbed actomyosin signaling. In addition, the underlying mechanism of HCM-induced hypercontractility, including R403Q HCM, has been explained in light of thick filament structural alterations and the tension generation (Alamo et al., 2017; Nag et al., 2017; Trivedi et al., 2018). Granted that, since with change in coefficient (Table 2) the XB cycling machinery could closely capture the HCM behavior (Figure 3), our model potentially attributes the HCM-induced disturbed interfilament signaling to a distortion-dependent forward transition from to state and the subsequent effect due to population of strongly bound XB states in the steady-state condition. This could imply that the HCM-induced disturbed actomyosin signaling might stem from a misregulated isomerization of myosin heads from pre-rotated to post-rotated force generating state affecting the strain induced in the myosin neck region. Analogically, the same process might explain how OM, BLEB, and Mava contribute to a new interfilament balance, restoring normal XB cycling tension generation. This hypothesis based on our model predictions could be insightful for accelerating the future drug development for sarcomere cardiomyopathies and mutation-specific HCM.

We have introduced a mechanistic solution to incorporate the dose-dependent effect of Blebbistatin and Omecamtiv mecarbil (Figure 4) consistent with experiments (Kampourakis et al., 2018) focusing on . Explicitly, the inverse Hill function trend observed in dose-dependent effects of OM and BLEB (Kampourakis et al., 2018) also reflects in the values of coefficients obtained for the studied concentrations (Figures 4B–E). Moreover, changes in F2, hf, and (Table 4) suggest that OM also favors the rapid detachment of XBs as another contributor to the disturbed actomyosin coupling involved in OM mechanism of action. Furthermore, the accurate simulation of shortened APD due to 1 µM OM combined with elevated basal ATPase rate values further signifies that inclusion of metabolite-sensitive transitions in the XB cycling of the CE cannot be ignored in precision medicine, especially when simulating the effects of drugs whose main mechanism of action impacts ATPases.

Limitations and future works

The developed mathematical model naturally has some limitations and potentials for advancements in the future studies. Firstly, as we use ODEs instead of computationally expensive PDEs, the cooperative spatial interactions between regulatory proteins and XB action have been approximated with a mean-field technique (Rice et al., 2008). Secondly, we assumed that 0.5 µM MAVA restores the elevated MgADP and Pi in HCM model variant to their basal values (listed in Table 2). Although this assumption is consistent with MAVA mechanism of action (Green et al., 2016), it is still a simplification. Thirdly, one fundamental limitation of hiPSC-CMs is that they rely on glycolytic metabolism, in contrast to the fatty acid-based metabolism of native human adult ventricular cardiomyocytes. A further source of energetic dissimilarities is the differences between surrounding medium in vitro vs. in vivo conditions. Given these and as our model does not include energy production process, we consider the metabolic differences out of the scope of this work. As detailed in vitro data on hiPSC-CM metabolism emerges, our modeling efforts serve as a solid basis for the next phase of cardiomyocyte models with energy production included.

The structural immaturity and special sarcomere alignment and performance in hiPSC-CMs and its effect on the HCM and drug-induced studies are important. The core of our HCM variant and MAVA calibration is based on hiPSC-CMs in vitro data obtained at day 30 post-differentiation, indicating cardiomyocytes maturation (Toepfer et al., 2020). In this work (Toepfer et al., 2020), to validate the hiPSC-CMs finding regarding in vivo data, the authors conducted parallel analyses on mouse HCM model and human HCM cardiac tissues. The comparison revealed that each HCM variant (HCM mice and hiPSC-CMs) caused hypercontractility with respect to its WT model (Figure 3 C&I in (Toepfer et al., 2020)). Furthermore, the prolonged contractile relaxation also was observed in mouse cardiomyocytes and hiPSC-CMs (Figures 3D,J in (Toepfer et al., 2020)). In addition, the dose-dependent decrease in hypercontractility and myosin population in DRX due to MAVA was observed in mouse cardiomyocytes and hiPSC-CMs (Figures 3E,K in (Toepfer et al., 2020)). All in all, Figure 3 in Toepfer et al. (2020) shows that HCM-induced and drug induced effects on hiPSC-CMs are consistent with the corresponding trend observed in mouse WT and R403Q cardiomyocytes variants. Although, our priority in calibration and validation of our results was using hiPSC-CMs in vitro data wherever available, all the abovementioned points imply that validation and calibration of HCM and drug-induced results with cell lines other than hiPSC-CMs are still valid and legitimate approaches.

Furthermore, as our computational model is 0D, it does not explicitly account for the disorganization of sarcomeres typical for hiPSC-CMs. This could be an interesting and valuable future direction for HCM computational studies as more experimental data becomes available.

Finally, the model can benefit from the inclusion of a metabolite-sensitive formulation of the intracellular SERCA pump, as a key ATP-dependent transporter. However, as the focus of this work was sarcomeric cardiomyopathies we have considered it out of scope here.

Conclusion

As cardiac precision medicine arises (Niederer et al., 2019; Paci et al., 2021; Forouzandehmehr et al., 2022), the demand for comprehensive computational models capable of performing high throughput pharmacological investigations heightens. This works proposes a novel metabolite-sensitive computational model of hiPSC-CMs electromechanics with demonstrated capacity to simulate sarcomeric cardiomyopathies and the compounds directly affecting the myosin dynamics considering the metabolic pathways. The mechanistic method offered for simulating the effects of HCM and drugs in this work lends insights upon the molecular interactions in contractile function and advance our pathophysiological understanding of the development of future therapeutics for HCM.

References

• 1

  AbrahamM. R.BottomleyP. A.DimaanoV. L.PinheiroA.SteinbergA.TraillT. A.et al (2013). Creatine kinase adenosine triphosphate and phosphocreatine energy supply in a single kindred of patients with hypertrophic cardiomyopathy. Am. J. Cardiol.112, 861–866. 10.1016/j.amjcard.2013.05.017

  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AlamoL.WareJ. S.PintoA.GillilanR. E.SeidmanJ. G.SeidmanC. E.et al (2017). Effects of myosin variants on interacting-heads motif explain distinct hypertrophic and dilated cardiomyopathy phenotypes. Elife6, e24634. 10.7554/eLife.24634

  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  AlsulamiK.MarstonS. (2020). Small molecules acting on myofilaments as treatments for heart and skeletal muscle diseases. Int. J. Mol. Sci.21, 9599. 10.3390/ijms21249599

  • CrossRef
  • Google Scholar
  • View reference in article

• 4

  AwindaP. O.BishawY.WatanabeM.GuglinM. A.CampbellK. S.TannerB. C. W. (2020). Effects of mavacamten on Ca 2+ sensitivity of contraction as sarcomere length varied in human myocardium. Br. J. Pharmacol.177, 5609–5621. 10.1111/bph.15271

  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  BakkehaugJ. P.KildalA. B.EngstadE. T.BoardmanN.NæsheimT.RønningL.et al (2015). Myosin activator omecamtiv mecarbil increases myocardial oxygen consumption and impairs cardiac efficiency mediated by resting myosin ATPase activity. Circ. Heart Fail.8, 766–775. 10.1161/CIRCHEARTFAILURE.114.002152

  • CrossRef
  • Google Scholar
  • View reference in article

• 6

  ForouzandehmehrM.GhoytasiI.ShamlooA.GhosiS. (2022). Particles in coronary circulation: A review on modelling for drug carrier design. Mat. Des.216, 110511. 10.1016/j.matdes.2022.110511

  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  ForouzandehmehrM.KoivumäkiJ. T.HyttinenJ.PaciM. (2021). A mathematical model of hiPSC cardiomyocytes electromechanics. Physiol. Rep.9, e15124. 10.14814/phy2.15124

  • CrossRef
  • Google Scholar
  • View reference in article

• 8

  GollapudiS. K.YuM.GanQ.-F.NagS. (2021). Synthetic thick filaments: A new avenue for better understanding the myosin super-relaxed state in healthy, diseased, and mavacamten-treated cardiac systems. J. Biol. Chem.296, 100114. 10.1074/jbc.RA120.016506

  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  GreenE. M.WakimotoH.AndersonR. L.EvanchikM. J.GorhamJ. M.HarrisonB. C.et al (2016). A small-molecule inhibitor of sarcomere contractility suppresses hypertrophic cardiomyopathy in mice. Science351, 617–621. 10.1126/science.aad3456

  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  GyimesiM.RauscherA. Á.SutharS. K.HamowK. Á.OraveczK.LőrinczI.et al (2021). Improved inhibitory and absorption, distribution, metabolism, excretion, and toxicology (ADMET) properties of blebbistatin derivatives indicate that blebbistatin scaffold is ideal for drug development targeting myosin-2. J. Pharmacol. Exp. Ther.376, 358–373. 10.1124/jpet.120.000167

  • CrossRef
  • Google Scholar
  • View reference in article

• 11

  HoC. Y.DayS. M.AshleyE. A.MichelsM.PereiraA. C.JacobyD.et al (2018). Genotype and lifetime burden of disease in hypertrophic cardiomyopathy: Insights from the sarcomeric human cardiomyopathy registry (SHaRe). Circulation138, 1387–1398. 10.1161/CIRCULATIONAHA.117.033200

  • CrossRef
  • Google Scholar
  • View reference in article

• 12

  IngwallJ. S. (2002). ATP and the heart. Boston, MA: Springer US. 10.1007/978-1-4615-1093-2

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  KampourakisT.ZhangX.SunY.-B.IrvingM. (2018). Omecamtiv mercabil and blebbistatin modulate cardiac contractility by perturbing the regulatory state of the myosin filament. J. Physiol.596, 31–46. 10.1113/JP275050

  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  KawasR. F.AndersonR. L.IngleS. R. B.SongY.SranA. S.RodriguezH. M. (2017). A small-molecule modulator of cardiac myosin acts on multiple stages of the myosin chemomechanical cycle. J. Biol. Chem.292, 16571–16577. 10.1074/jbc.M117.776815

  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  MalhotraA.SharmaS. (2017). Ischemic heart disease in women: Not about religion. Eur. Cardiol.12, 8–9. 10.15420/ecr.2017.12.1.GE1

  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  MargaraF.RodriguezB.ToepferC. N.Bueno-OrovioA. (2021). Mavacamten efficacy in mutation-specific hypertrophic cardiomyopathy: an in silico approach to inform precision medicine. Comput. Cardiol.2021, 1–4. 10.23919/CinC53138.2021.9662736

  • CrossRef
  • Google Scholar
  • View reference in article

• 17

  NagS.SommeseR. F.UjfalusiZ.CombsA.LangerS.SuttonS.et al (2015). Contractility parameters of human β-cardiac myosin with the hypertrophic cardiomyopathy mutation R403Q show loss of motor function. Sci. Adv.1, e1500511. 10.1126/sciadv.1500511

  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  NagS.V TrivediD.SarkarS. S.AdhikariA. S.SunithaM. S.SuttonS.et al (2017). The myosin mesa and the basis of hypercontractility caused by hypertrophic cardiomyopathy mutations. Nat. Struct. Mol. Biol.24, 525–533. 10.1038/nsmb.3408

  • CrossRef
  • Google Scholar
  • View reference in article

• 19

  NiedererS. A.LumensJ.TrayanovaN. A. (2019). Computational models in cardiology. Nat. Rev. Cardiol.16, 100–111. 10.1038/s41569-018-0104-y

  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  PaciM.KoivumäkiJ. T.LuH. R.GallacherD. J.PassiniE.RodriguezB. (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, 604713. 10.3389/fphar.2021.604713

  • CrossRef
  • Google Scholar
  • View reference in article

• 21

  PaciM.PassiniE.KlimasA.SeveriS.HyttinenJ.RodriguezB.et al (2020). All-optical electrophysiology refines populations of in silico human iPSC-CMs for drug evaluation. Biophys. J.118, 2596–2611. 10.1016/j.bpj.2020.03.018

  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  PaciM.PölönenR. P.CoriD.PenttinenK.Aalto-SetäläK.SeveriS.et al (2018). Automatic optimization of an in silico model of human iPSC derived cardiomyocytes recapitulating calcium handling abnormalities. Front. Physiol.9, 709. 10.3389/fphys.2018.00709

  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  PionerJ. M.GuanX.KlaimanJ. M.RaccaA. W.PabonL.MuskheliV.et al (2020). Absence of full-length dystrophin impairs normal maturation and contraction of cardiomyocytes derived from human-induced pluripotent stem cells. Cardiovasc. Res.116, 368–382. 10.1093/cvr/cvz109

  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  QuY.GaoB.ArimuraZ.FangM.VargasH. M. (2021). Comprehensive in vitro pro‐arrhythmic assays demonstrate that omecamtiv mecarbil has low pro‐arrhythmic risk. Clin. Transl. Sci.14, 1600–1610. 10.1111/cts.13039

  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  RahmanM. A.UšajM.RassierD. E.MånssonA. (2018). Blebbistatin effects expose hidden secrets in the force-generating cycle of actin and myosin. Biophys. J.115, 386–397. 10.1016/j.bpj.2018.05.037

  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  RiceJ. J.WangF.BersD. M.De TombeP. P. (2008). Approximate model of cooperative activation and crossbridge cycling in cardiac muscle using ordinary differential equations. Biophys. J.95, 2368–2390. 10.1529/biophysj.107.119487

  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  RohdeJ. A.RoopnarineO.ThomasD. D.MurettaJ. M. (2018). Mavacamten stabilizes an autoinhibited state of two-headed cardiac myosin. Proc. Natl. Acad. Sci. U. S. A.115, E7486–E7494. 10.1073/pnas.1720342115

  • CrossRef
  • Google Scholar

• 28

  RuanJ. L.TullochN. L.RazumovaM. V.SaigetM.MuskheliV.PabonL.et al (2016). Mechanical stress conditioning and electrical stimulation promote contractility and force maturation of induced pluripotent stem cell-derived human cardiac tissue. Circulation134, 1557–1567. 10.1161/CIRCULATIONAHA.114.014998

  • CrossRef
  • Google Scholar

• 29

  SantiniL.PalandriC.NedianiC.CerbaiE.CoppiniR. (2020). Modelling genetic diseases for drug development: Hypertrophic cardiomyopathy. Pharmacol. Res.160, 105176. 10.1016/j.phrs.2020.105176

  • CrossRef
  • Google Scholar

• 30

  SarkarS. S.TrivediD. V.MorckM. M.AdhikariA. S.PashaS. N.RuppelK. 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, eaax0069. 10.1126/sciadv.aax0069

  • CrossRef
  • Google Scholar

• 31

  SchmidM.ToepferC. N. (2021). Cardiac myosin super relaxation (SRX): a perspective on fundamental biology, human disease and therapeutics. Biol. Open10, bio057646. 10.1242/bio.057646

  • CrossRef
  • Google Scholar

• 32

  SemsarianC.InglesJ.MaronM. S.MaronB. J. (2015). New perspectives on the prevalence of hypertrophic cardiomyopathy. J. Am. Coll. Cardiol.65, 1249–1254. 10.1016/j.jacc.2015.01.019

  • CrossRef
  • Google Scholar

• 33

  SequeiraV.BerteroE.MaackC. (2019). Energetic drain driving hypertrophic cardiomyopathy. FEBS Lett.593, 1616–1626. 10.1002/1873-3468.13496

  • CrossRef
  • Google Scholar

• 34

  SewananL. R.SchwanJ.KlugerJ.ParkJ.JacobyD. L.QyangY.et al (2019). Extracellular matrix from hypertrophic myocardium provokes impaired twitch dynamics in healthy cardiomyocytes. JACC. Basic Transl. Sci.4, 495–505. 10.1016/j.jacbts.2019.03.004

  • CrossRef
  • Google Scholar

• 35

  SpindlerM.SaupeK. W.ChristeM. E.SweeneyH. L.SeidmanC. E.SeidmanJ. G.et al (1998). Diastolic dysfunction and altered energetics in the alphaMHC403/+ mouse model of familial hypertrophic cardiomyopathy. J. Clin. Invest.101, 1775–1783. 10.1172/JCI1940

  • CrossRef
  • Google Scholar

• 36

  SpudichJ. A. (2014). Hypertrophic and dilated cardiomyopathy: four decades of basic research on muscle lead to potential therapeutic approaches to these devastating genetic diseases. Biophys. J.106, 1236–1249. 10.1016/j.bpj.2014.02.011

  • CrossRef
  • Google Scholar

• 37

  SzentandrassyN.HorvathB.VacziK.KistamasK.MasudaL.MagyarJ.et al (2016). Dose-dependent electrophysiological effects of the myosin activator omecamtiv mecarbil in canine ventricular cardiomyocytes. J. Physiol. Pharmacol.67, 483–489.

  • Google Scholar

• 38

  ToepferC. N.GarfinkelA. C.VenturiniG.WakimotoH.RepettiG.AlamoL.et al (2020). Myosin sequestration regulates sarcomere function, cardiomyocyte energetics, and metabolism, informing the pathogenesis of hypertrophic cardiomyopathy. Circulation141, 828–842. 10.1161/CIRCULATIONAHA.119.042339

  • CrossRef
  • Google Scholar

• 39

  ToepferC. N.WakimotoH.GarfinkelA. C.McDonoughB.LiaoD.JiangJ.et al (2019). Hypertrophic cardiomyopathy mutations in MYBPC3 dysregulate myosin. Sci. Transl. Med.11, eaat1199. 10.1126/scitranslmed.aat1199

  • CrossRef
  • Google Scholar

• 40

  TranK.HanJ.-C.CrampinE. J.TabernerA. J.LoiselleD. S. (2017). Experimental and modelling evidence of shortening heat in cardiac muscle. J. Physiol.595, 6313–6326. 10.1113/JP274680

  • CrossRef
  • Google Scholar

• 41

  TranK.SmithN. P.LoiselleD. S.CrampinE. J. (2010). A metabolite-sensitive, thermodynamically constrained model of cardiac cross-bridge cycling: Implications for force development during ischemia. Biophys. J.98, 267–276. 10.1016/j.bpj.2009.10.011

  • CrossRef
  • Google Scholar

• 42

  TrivediD. V.AdhikariA. S.SarkarS. S.RuppelK. M.SpudichJ. A. (2018). Hypertrophic cardiomyopathy and the myosin mesa: viewing an old disease in a new light. Biophys. Rev.10, 27–48. 10.1007/s12551-017-0274-6

  • CrossRef
  • Google Scholar

• 43

  TsukamotoO. (2019). Direct sarcomere modulators are promising new treatments for cardiomyopathies. Int. J. Mol. Sci.21, 226. 10.3390/ijms21010226

  • CrossRef
  • Google Scholar

• 44

  UtterM. S.RybaD. M.LiB. H.WolskaB. M.SolaroR. J. (2015). Omecamtiv mecarbil, a cardiac myosin activator, increases Ca2+ sensitivity in myofilaments with a dilated cardiomyopathy mutant tropomyosin E54K. J. Cardiovasc. Pharmacol.66, 347–353. 10.1097/FJC.0000000000000286

  • CrossRef
  • Google Scholar

• 45

  van der VeldenJ.TocchettiC. G.VarricchiG.BiancoA.SequeiraV.Hilfiker-KleinerD.et al (2018). Metabolic changes in hypertrophic cardiomyopathies: scientific update from the working group of myocardial function of the European society of cardiology. Cardiovasc. Res.114, 1273–1280. 10.1093/cvr/cvy147

  • CrossRef
  • Google Scholar

• 46

  WangL.KimK.ParikhS.CadarA. G.BersellK. R.HeH.et al (2018). Hypertrophic cardiomyopathy-linked mutation in troponin T causes myofibrillar disarray and pro-arrhythmic action potential changes in human iPSC cardiomyocytes. J. Mol. Cell. Cardiol.114, 320–327. 10.1016/j.yjmcc.2017.12.002

  • CrossRef
  • Google Scholar