- Faculty of Medicine and Health Technology, Tampere University, Tampere, Finland
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 Ca2+ 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-Ca2+ 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 Ca2+ 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-Ca2+ 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 IC50 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 Ca2+ sensitivity, ATPase dynamics, and dose-dependent effect of MAVA, BLEB, and OM on the tension-Ca2+ 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 Ca2+ sensitivity, contractility, and energetics of the hiPSC-CMs.
FIGURE 1. Schematics of the interfilament coupling in cardiac force generation (A) the modelled crossbridge (XB) cycling used in hiMCE model (B), and the schematic of hiPSC-CM cell main functional components (C). DRX: Disturbed relax state. SRX: Super relaxed state. T: troponin, TCa: Ca2+ bound troponin, NXB: non-permissive state preventing XB formation, PXB: permissive state of XB formation, XBpreR: strongly bound XB before isomerised rotation, XBpstR: XB in strongly bound post isomerised rotation state, AM1 and AM2 are strongly-bound rapid equilibrium substates contributing equally to the force generation and we assumed MgADP binds to AM1 (Tran et al., 2010).
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 Ca2+ 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 (XBpstR) 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 Ca2+ activation is divided into four states including a non-permissive (NXB), a permissive (PXB), a pre power stroke state (XBpreR), and a post power stroke state (XBpstR) (Figure 1) (Tran et al., 2017). XBpreR and XBpstR 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 NXB state and when activated by Ca2+ XBs enter PXB 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.
TABLE 1. The values of the contractile element parameters for hiPSC-CM-CE (Forouzandehmehr et al., 2021) and hiMCE. F1 and F2 represent DRX:SRX/(DRX:SRX)control ratio in the XB cycling. Kon, Knp and Kpn are rate constants for Ca2+ binding to troponin, the forward and backward transition rates between Pxb and Nxb states, respectively (Rice et al., 2008). Nperm and perm50 denote the Hill coefficient and the half-activation constant, respectively, describing the nonlinearity of the cooperativity in Ca2+ activation of XBs (Rice et al., 2008). KoffL and KoffH represent the rate constants affecting Ca2+ unbinding from low and high affinity sites on troponin, respectively (Rice et al., 2008). M denotes the mass term in the model of sarcomere by Rice et al. (Rice et al., 2008). Kxb represents the tension scaler detailed in (Forouzandehmehr et al., 2021). Xbmodsp is a species-dpendent XB cycling rate scaler (Rice et al., 2008). Hf denotes the rate constant in the forward transition between XBpreR and XBpstR (Rice et al., 2008).
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:
Where
Again,
The value for
Furthermore,
Where
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).
TABLE 2. The parameter used in the HCM model variant. Pi_ref denotes the reference value for inorganic phosphate (Pi) in the simulations. Ap2 is a variable influencing detachment of crossbridges. F1 and F2 are coefficients affecting pre-rotational states in XB cycling as also used in (Margara et al., 2021).
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 XBpreR and PXB 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 XBpreR 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 Ca2+ activation and binding process (Awinda et al., 2020), Pi (Alsulami and Marston, 2020), and ATPase activity (Green et al., 2016). Specifically, changes in Kon and nperm 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,
TABLE 3. Modifications to the model parameters to simulate the effect of Mavacamten. A, B, C, D, and E are coefficients in Eqs 15–19 affecting PXB-XBPreR regulation, Pi-dependent transition in PXB-XBPreR, XBpreR-XBpostR regulation, Proton-dependent transition in XBpostR-XBpreR, and MgATP-dependent transition from XBpostR-PXB, respectively. BL is the baseline value given in Table 1. Default values of A-E, ap1 and ap3 coefs., F1 and F2 are equal to 1.
To indicate, A modulates PXB to XBpreR transition, B influences Pi-dependent XBpreR to PXB transition, C takes effect on XBpreR to XBpstR transition, D affects proton-dependent XBpstR to XBpreR transition, and E controls MgATP-dependent transition from XBpstR to PXB 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
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,
TABLE 4. The modifications to the parameters to simulate the effects of BLEB and OM. BL: Baseline values of hiMCE model given in Table 1.
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,
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.
TABLE 5. The experimental data used for calibration of the model and validation of the simulated results.
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, Ca2+ 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).
TABLE 6. Action potential (AP), Ca2+ transients (CaT), and active tension (AT) calculated biomarkers in spontaneous condition and their comparison with Paci2020 and hiPSC-CM-CE model (i.e. no metabolite-sensitive CE) and the experimental values (Paci et al., 2018; Paci et al., 2020). APA: AP amplitude, MDP: maximum diastolic potential, CL: cycle length, dV/dt max: maximum upstroke velocity, APD10 and APD30 and APD90: AP duration at 10, 30, 90% of repolarization, respectively, AP Tri: AP triangulation index. The simulated biomarkers of CaT are DURATION: duration of the transient, tRise10, peak: time to peak, tRise10, 50 and tRise10, 90: rise time from 10 to 50% and 90% of maximum threshold, respectively, and tDecay90,10: decay time from 90 to 10%. AT: Active tension, RT50: time from peak contraction to 50% of relaxation, %FS: percent of fractional shortening. The experimental ranges for contraction biomarkers are from (Forouzandehmehr et al., 2021). The third column is taken directly from the original Paci2020 publication (Paci et al., 2020).
FIGURE 2. Standard results of the model in spontaneous beating: Action Potentials (A), Ca2+ Transients (B), Active Tensions (C), ATPase rate (D), Flux of Ca2+ towards the contractile element (E), Fractional cell shortening at 1 Hz pacing (F). Cited works: (Ruan et al., 2016; 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. Simulated action potential (A), Calcium transients (B), active tensions (C), ATPase rate (D), Flux of Ca2+ towards the myofilament (E), and fractional cell shortenings (F) in R403Q hypertrophic cardiomyopathy and Mavacamten modes (All simulations were done at 1 Hz pacing). The percents of prolonged tension relaxation in R403Q mode (C), the reduction in tension relaxation due to MAVA (C), the reduction in fractional shortening in R403Q mode due to MAVA (F), and the reduction in ATPase rate due to MAVA (D) agree with the experimental data (Toepfer et al., 2020; 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 pCa50 in the tension-Ca2+ 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 Ca2+ 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,
TABLE 7. Experimental data (Kampourakis et al., 2018) and hiMCE results due to the effect of 1 µM OM and 5 µM BLEB. pCa50 represents the -log of the Ca2+ concentration associated with 50% of maximum tension. cTnC-E: data from cardiac troponin C (cTnC) E-helix obtained by a rhodamine probe. cRLC-E: Data from a probe connected to the myosin regulatory light chain (RLC).
FIGURE 4. Calcium-tension relationships in control and drug-induced modes in isometric condition (A), ap2 coefficients found for different BLEB and OM concentrations (B,C), and dose-dependent tension-Ca2+ relationships of BLEB (D) and OM (E) in isometric conditions. OM: Omecamtive mecarbil. BLEB: Blebbistatin. Experimental data from (Kampourakis et al., 2018). T0: 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 APD90 values (Figure 5A). This translates to 17 ms reduction in APD90 (502–485 ms) which is consistent with the order of APD90 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. Predicted effect of 1 µM OM by hiMCE on action potentials (A), Ca2+ flux towards the myofilament (B), L-type Ca2+ current (C), Na+/Ca2+ exchanger INCX (D) Ca2+ transients (E), sarcolemmal Ca2+ pump current (F), Na+/K+ pump (G), and ATPase rates (H). The change of APD has been considered calculating APD90 in agreement with experimental data reported in (Szentandrassy et al., 2016). OM: Omecamtiv mecarbil.
Our model predicts a 23% increase in the amplitude of the Ca2+ flux towards the myofilament, JCB (Figure 5B). The subsequent accumulation of intracellular Ca2+ 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 ICaL (Figure 5C) and only very slightly increased IpCa (Figure 5F). Whereas there is a more substantial 17% increase in the amplitude of the INCX (Figure 5D). The enhanced reverse mode of INCX causes accumulation of intracellular Na+ (7.04 vs. 7.43 mM) that promotes a stronger repolarizing INaK (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 (PXB) and the strongly bound XBs before isomerized rotation (XBpreR) state. Also, we modulated Ca2+ binding and sensitivity of the CE, and MgATP-dependent transitions between PXB and XBs in strongly bound post isomerized rotation state (XBpostR) 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 XBpreR and XBpostR 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 Ca2+ 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 koff constants which are the rate constants affecting Ca2+ 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 Ca2+ 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 Kon, Knp and Kpn rate constants for Ca2+ binding to troponin, forward and backward transition rates between Pxb and Nxb states, respectively (Table 4). Further, the constant rates directly affecting Ca2+ bound troponin thin filaments regulation induced by Ca2+ bound troponin (Kon, nperm, koffL, koffH, kpn, knp, and koffmode) 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 Ca2+ 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 Ca2+ 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 Ca2+ 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
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
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.
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
MF, MP, JK and JH contributed to conception and design of the study. MF performed the simulations and wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.
Funding
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. JK was supported by Academy of Finland Centre of Excellence in Body-on-Chip Research, Pirkanmaa regional fund of the Finnish Cultural Foundation (grant number 50201322), and Finnish Foundation for Cardiovascular Research (grant number 200101).
Acknowledgments
The authors wish to thank Samuel Wall for very helpful comments during the preparation of this manuscript.
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.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2022.1010786/full#supplementary-material
References
Abraham M. R., Bottomley P. A., Dimaano V. L., Pinheiro A., Steinberg A., Traill T. 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. doi:10.1016/j.amjcard.2013.05.017
Alamo L., Ware J. S., Pinto A., Gillilan R. E., Seidman J. G., Seidman C. E., et al. (2017). Effects of myosin variants on interacting-heads motif explain distinct hypertrophic and dilated cardiomyopathy phenotypes. Elife 6, e24634. doi:10.7554/eLife.24634
Alsulami K., Marston S. (2020). Small molecules acting on myofilaments as treatments for heart and skeletal muscle diseases. Int. J. Mol. Sci. 21, 9599. doi:10.3390/ijms21249599
Awinda P. O., Bishaw Y., Watanabe M., Guglin M. A., Campbell K. S., Tanner B. 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. doi:10.1111/bph.15271
Bakkehaug J. P., Kildal A. B., Engstad E. T., Boardman N., Næsheim T., Rønning L., 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. doi:10.1161/CIRCHEARTFAILURE.114.002152
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
Forouzandehmehr M., Koivumäki J. T., Hyttinen J., Paci M. (2021). A mathematical model of hiPSC cardiomyocytes electromechanics. Physiol. Rep. 9, e15124. doi:10.14814/phy2.15124
Gollapudi S. K., Yu M., Gan Q.-F., Nag S. (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. doi:10.1074/jbc.RA120.016506
Green E. M., Wakimoto H., Anderson R. L., Evanchik M. J., Gorham J. M., Harrison B. C., et al. (2016). A small-molecule inhibitor of sarcomere contractility suppresses hypertrophic cardiomyopathy in mice. Science 351, 617–621. doi:10.1126/science.aad3456
Gyimesi M., Rauscher A. Á., Suthar S. K., Hamow K. Á., Oravecz K., Lőrincz I., 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. doi:10.1124/jpet.120.000167
Ho C. Y., Day S. M., Ashley E. A., Michels M., Pereira A. C., Jacoby D., et al. (2018). Genotype and lifetime burden of disease in hypertrophic cardiomyopathy: Insights from the sarcomeric human cardiomyopathy registry (SHaRe). Circulation 138, 1387–1398. doi:10.1161/CIRCULATIONAHA.117.033200
Kampourakis T., Zhang X., Sun Y.-B., Irving M. (2018). Omecamtiv mercabil and blebbistatin modulate cardiac contractility by perturbing the regulatory state of the myosin filament. J. Physiol. 596, 31–46. doi:10.1113/JP275050
Kawas R. F., Anderson R. L., Ingle S. R. B., Song Y., Sran A. S., Rodriguez H. M. (2017). A small-molecule modulator of cardiac myosin acts on multiple stages of the myosin chemomechanical cycle. J. Biol. Chem. 292, 16571–16577. doi:10.1074/jbc.M117.776815
Malhotra A., Sharma S. (2017). Ischemic heart disease in women: Not about religion. Eur. Cardiol. 12, 8–9. doi:10.15420/ecr.2017.12.1.GE1
Margara F., Rodriguez B., Toepfer C. N., Bueno-Orovio A. (2021). Mavacamten efficacy in mutation-specific hypertrophic cardiomyopathy: an in silico approach to inform precision medicine. Comput. Cardiol. 2021, 1–4. doi:10.23919/CinC53138.2021.9662736
Nag S., Sommese R. F., Ujfalusi Z., Combs A., Langer S., Sutton S., et al. (2015). Contractility parameters of human β-cardiac myosin with the hypertrophic cardiomyopathy mutation R403Q show loss of motor function. Sci. Adv. 1, e1500511. doi:10.1126/sciadv.1500511
Nag S., V Trivedi D., Sarkar S. S., Adhikari A. S., Sunitha M. S., Sutton S., et al. (2017). The myosin mesa and the basis of hypercontractility caused by hypertrophic cardiomyopathy mutations. Nat. Struct. Mol. Biol. 24, 525–533. doi:10.1038/nsmb.3408
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
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, 604713. doi:10.3389/fphar.2021.604713
Paci M., Passini E., Klimas A., Severi S., Hyttinen J., Rodriguez B., et al. (2020). All-optical electrophysiology refines populations of in silico human iPSC-CMs for drug evaluation. Biophys. J. 118, 2596–2611. doi:10.1016/j.bpj.2020.03.018
Paci M., Pölönen R. P., Cori D., Penttinen K., Aalto-Setälä K., Severi S., et al. (2018). Automatic optimization of an in silico model of human iPSC derived cardiomyocytes recapitulating calcium handling abnormalities. Front. Physiol. 9, 709. doi:10.3389/fphys.2018.00709
Pioner J. M., Guan X., Klaiman J. M., Racca A. W., Pabon L., Muskheli V., 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. doi:10.1093/cvr/cvz109
Qu Y., Gao B., Arimura Z., Fang M., Vargas H. M. (2021). Comprehensive in vitro pro‐arrhythmic assays demonstrate that omecamtiv mecarbil has low pro‐arrhythmic risk. Clin. Transl. Sci. 14, 1600–1610. doi:10.1111/cts.13039
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
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
Rohde J. A., Roopnarine O., Thomas D. D., Muretta J. M. (2018). Mavacamten stabilizes an autoinhibited state of two-headed cardiac myosin. Proc. Natl. Acad. Sci. U. S. A. 115, E7486–E7494. doi:10.1073/pnas.1720342115
Ruan J. L., Tulloch N. L., Razumova M. V., Saiget M., Muskheli V., Pabon L., et al. (2016). Mechanical stress conditioning and electrical stimulation promote contractility and force maturation of induced pluripotent stem cell-derived human cardiac tissue. Circulation 134, 1557–1567. doi:10.1161/CIRCULATIONAHA.114.014998
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
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, eaax0069. doi:10.1126/sciadv.aax0069
Schmid M., Toepfer C. N. (2021). Cardiac myosin super relaxation (SRX): a perspective on fundamental biology, human disease and therapeutics. Biol. Open 10, bio057646. doi:10.1242/bio.057646
Semsarian C., Ingles J., Maron M. S., Maron B. J. (2015). New perspectives on the prevalence of hypertrophic cardiomyopathy. J. Am. Coll. Cardiol. 65, 1249–1254. doi:10.1016/j.jacc.2015.01.019
Sequeira V., Bertero E., Maack C. (2019). Energetic drain driving hypertrophic cardiomyopathy. FEBS Lett. 593, 1616–1626. doi:10.1002/1873-3468.13496
Sewanan L. R., Schwan J., Kluger J., Park J., Jacoby D. L., Qyang Y., et al. (2019). Extracellular matrix from hypertrophic myocardium provokes impaired twitch dynamics in healthy cardiomyocytes. JACC. Basic Transl. Sci. 4, 495–505. doi:10.1016/j.jacbts.2019.03.004
Spindler M., Saupe K. W., Christe M. E., Sweeney H. L., Seidman C. E., Seidman J. G., et al. (1998). Diastolic dysfunction and altered energetics in the alphaMHC403/+ mouse model of familial hypertrophic cardiomyopathy. J. Clin. Invest. 101, 1775–1783. doi:10.1172/JCI1940
Spudich J. 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. doi:10.1016/j.bpj.2014.02.011
Szentandrassy N., Horvath B., Vaczi K., Kistamas K., Masuda L., Magyar J., et al. (2016). Dose-dependent electrophysiological effects of the myosin activator omecamtiv mecarbil in canine ventricular cardiomyocytes. J. Physiol. Pharmacol. 67, 483–489.
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
Toepfer C. N., Wakimoto H., Garfinkel A. C., McDonough B., Liao D., Jiang J., et al. (2019). Hypertrophic cardiomyopathy mutations in MYBPC3 dysregulate myosin. Sci. Transl. Med. 11, eaat1199. doi:10.1126/scitranslmed.aat1199
Tran K., Han J.-C., Crampin E. J., Taberner A. J., Loiselle D. S. (2017). Experimental and modelling evidence of shortening heat in cardiac muscle. J. Physiol. 595, 6313–6326. doi:10.1113/JP274680
Tran K., Smith N. P., Loiselle D. S., Crampin E. J. (2010). A metabolite-sensitive, thermodynamically constrained model of cardiac cross-bridge cycling: Implications for force development during ischemia. Biophys. J. 98, 267–276. doi:10.1016/j.bpj.2009.10.011
Trivedi D. V., Adhikari A. S., Sarkar S. S., Ruppel K. M., Spudich J. A. (2018). Hypertrophic cardiomyopathy and the myosin mesa: viewing an old disease in a new light. Biophys. Rev. 10, 27–48. doi:10.1007/s12551-017-0274-6
Tsukamoto O. (2019). Direct sarcomere modulators are promising new treatments for cardiomyopathies. Int. J. Mol. Sci. 21, 226. doi:10.3390/ijms21010226
Utter M. S., Ryba D. M., Li B. H., Wolska B. M., Solaro R. 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. doi:10.1097/FJC.0000000000000286
van der Velden J., Tocchetti C. G., Varricchi G., Bianco A., Sequeira V., Hilfiker-Kleiner D., 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. doi:10.1093/cvr/cvy147
Wang L., Kim K., Parikh S., Cadar A. G., Bersell K. R., He H., 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. doi:10.1016/j.yjmcc.2017.12.002
Keywords: in silico modeling, human stem cell-derived cardiomyocyte, action potential, immature cardiomyocytes, cardiac metabolism, hypertrophic cardiomyopathy, pharmacology
Citation: Forouzandehmehr M, Paci M, Koivumäki JT and Hyttinen J (2022) Altered contractility in mutation-specific hypertrophic cardiomyopathy: A mechano-energetic in silico study with pharmacological insights. Front. Physiol. 13:1010786. doi: 10.3389/fphys.2022.1010786
Received: 03 August 2022; Accepted: 14 October 2022;
Published: 31 October 2022.
Edited by:
Simone Scacchi, University of Milan, ItalyReviewed by:
Serdar Göktepe, Middle East Technical University, TurkeyGuido Caluori, INSERM Institut de Rythmologie et Modélisation Cardiaque (IHU-Liryc), France
Copyright © 2022 Forouzandehmehr, Paci, Koivumäki and Hyttinen. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Mohamadamin Forouzandehmehr, mohamadamin.forouzandehmehr@tuni.fi