- Original Paper
A new myofilament contraction model with ATP consumption for ventricular cell model
The Journal of Physiological Sciences volume 68, pages 541–554 (2018)
A new contraction model of cardiac muscle was developed by combining previously described biochemical and biophysical models. The biochemical component of the new contraction model represents events in the presence of Ca2+–crossbridge attachment and power stroke following inorganic phosphate release, detachment evoked by the replacement of ADP by ATP, ATP hydrolysis, and recovery stroke. The biophysical component focuses on Ca2+ activation and force (F b) development assuming an equivalent crossbridge. The new model faithfully incorporates the major characteristics of the biochemical and biophysical models, such as F b activation by transient Ca2+ ([Ca2+]–F b), [Ca2+]–ATP hydrolysis relations, sarcomere length–F b, and F b recovery after jumps in length under the isometric mode and upon sarcomere shortening after a rapid release of mechanical load under the isotonic mode together with the load–velocity relationship. ATP consumption was obtained for all responses. When incorporated in a ventricular cell model, the contraction model was found to share approximately 60% of the total ATP usage in the cell model.
Activation of muscle contraction results in a cycle of crossbridge (CB) attachment–detachment between the two major sliding filaments, actin and myosin. During this cycle, the CB converts the biochemical energy to contractile mechanical force by hydrolyzing ATP catalyzed by actomyosin (AM) ATPase. Since this is the major component of ATP consumption in the cardiac muscle , a detailed mathematical model of CB dynamics and associated ATPase activity is a prerequisite for analyzing the energetics of cardiac muscle contraction under physiological conditions.
A variety of whole-cell models have been developed to analyze mechanisms underlying the physiological regulation of cardiac contraction. Recently we developed a human ventricular cell model that includes both a detailed calcium ion (Ca2+)-induced Ca2+ release model and a contraction model in addition to membrane excitation supported by the regulation of intracellular ion concentrations [2, 3]. However, to our knowledge, no cell model has been reported to date which calculates ATP consumption during electrical excitation and contraction of the muscle based on the detailed reaction cascade of ATP hydrolysis by the S1 segment of myosin.
Recent studies have clarified many details of the biochemical and structural states of the CB cycle, and 11 predominant reaction states of AM ATPase have been identified in skeletal muscle. Lombardi and Piazzesi , Månsson , and Månsson et al.  reduced the scheme to four to six essential lumped states. These myofilament models are based on the Huxley hypothesis of independent CB behavior and calculate the state transition of CB based on the hypothetical Gibb’s free energy profile as a function of relative distance between myosin head and the nearest actin binding site for each CB (spatially explicit model). This type of model has been continuously improved since the original model was proposed (see review by Månsson et al. ). However, the model is still difficult to use in the whole-cell model of excitation–contraction coupling because the model is based on CB kinetics in the continuous presence of a steady Ca2+ concentration [Ca2+] during tetanus contraction of skeletal muscle. Moreover, the profile of free energy of various states of CB is still largely variable between the different models that have been proposed [4, 5, 7], especially when the critical parameter of sliding velocity of the myofilament is included in the kinetics.
The objective of the study reported here is to develop a new hybrid model which combines the reaction steps of ATP hydrolysis by the S1 segment (‘Månsson model’) with the biophysical contraction models that describe the Ca2+-mediated activation of the myofilaments as well as the CB kinetics using mean-field approximation [8,9,10,11,12,13,14] instead of spatial explicit approaches [15, 16]. The Negroni and Lascano model (‘NL model’ [13, 14]) assumes a relatively simple equivalent CB (eCB), the distortion of which is represented with a linear spring and its sliding along the actin filament calculated by assuming a movable viscosity head. These approximations seem to be roughly in common with the Rice model . The NL model has also been used to develop various types of ventricular models. Thus, we refer to the NL model in the present study because of its utility in developing the hybrid model. The ‘NL model’ successfully reconstructs most of the key mechanical events of the cardiac muscle described to date in biophysical studies, such as the time course of developing the tension evoked by the membrane excitation, as well as the tensions evoked by applying step changes in the isometric and isotonic contraction modes at various [Ca2+]. The model also successfully reconstructs the relationships of [Ca2+] and force, of force and half-sarcomere length, and of force and velocity of shortening [17, 18].
The development of such a hybrid model may critically depend on reconstruction of the rapid time course of responses to square pulse application during the isometric or isotonic contraction modes, since such experiments directly measure the rates of CB movement as well as the state transition of the myofilament reactions in situ. The original rate constants, which were theoretically elaborated in the Månsson model, were simplified in reference to the kinetics of eCB. Thus, in the present study, we readjusted the rate of state transitions to obtain a new contraction model. When the hybrid contraction model was incorporated into the human ventricular model , ATP consumption by AM ATPase was approximately 60% of the total usage of the cell model. ATP consumption proportional to the pressure–volume trajectory area was obtained when the myofilament model was used in the ‘Laplace heart’ in the multi-scale hemodynamic model of cardiovascular system .
Materials and methods
The model of CB dynamics in the NL model is characterized by the introduction of the eCB, which largely facilitates the calculation of the developed force of contraction in a continuous manner. Even though the simplification was achieved without truly developing a specific algorithm for averaging the behavior of all CBs, the magnitude as well as the time-course of the reconstructed F b agrees well with the experimental findings . This eCB assumption was adopted in the present model as in the original study. Details on kinetic equations in the NL model are described in the Electronic Supplementary Materials (ESM), with the exception of those described in the following subsections. For details on the rate constants defined in the Månsson model, the reader is referred to the original publications [4,5,6].
Model reaction scheme
The lumped four-state biochemical model developed by Månsson  well represents the four essential steps of ATP hydrolysis. In the new state–transition scheme shown in Fig. 1, these four states are represented by symbols of cAM DPw, cAM Ds, cAM T, and cM DP (cM T), respectively, where A, M, T, D, P, and c refer to actin, myosin, ATP, ADP, inorganic phosphate (Pi), and Ca2+, respectively, and the lowercase letters w and s indicate weakly and strongly bound conformations of CB, respectively. The correspondence of each state with the original A 0 , A 1 , A 2 and A 3 states defined in the Månsson model is indicated below each state in the scheme shown in Fig. 1. A schematic composition of each state is indicated in the right panel with supplemental troponin–Ca2+ binding. Note that calmodulin is aligned on the thin filament, and the Ca2+ binding sites are occupied in all cM DP, cAM DPw, cAM Ds, cAM T, and cM DP (cM T), in contrast to the three unoccupied states, AM Ds, AM T, M DP (M T) aligned on the left relaxation pathway and associated with the same rate constants, Z d and Z h, respectively, for simplicity. The angle of the lever arm in connection to the S1 segment shows the conformation of the myosin lever arm, i.e., the power stroke at cAM Ds and the recovery stroke at cM DP and M DP. Thus, the Ca2+-binding steps of Q 1–Q 2 and Q 10–Q 11 as well as the CB attachment steps Q 3–Q 4 are adopted from the NL model, while the other steps are all based on the Månsson model.
Rate constants for these four Ca2+-bound states of ATP hydrolysis were adjusted in the present study to approximate the overall rate of ATP hydrolysis in situ at physiological temperature [1, 20], with reference to the theoretical ones determined as a function of ΔG of each conformation in the original Månsson model. The NL model is also well adapted to the kinetics of the ventricular myocardium at 37 °C. The time courses of the Ca2+ transition states followed by CB association and dissociation are quite similar to those observed in experiments ; namely, at systolic free intracellular Ca2+ ([Ca2+]i) of approximately 0.5–1 μM at physiological temperature, the increase in [Ca2+]i and the Ca2+-induced troponin-switch, which corresponds to the M DP → cM DP step, takes places within approximately 20 ms [22, 23].
Rate function for state transitions in the NL model
In the original NL model , the state transition from M DP to cM DP represents lumped steps of the Ca2+ binding to troponin and the subsequent conformation change in the tropomyosin complex. Therefore, in the new NL model this step was conserved as in the original model. In the original model, a hypothetical troponin system (TS) was assumed to consist of three sequential troponin–tropomyosin regulatory units, assuming a nearest neighbor influence in the Ca2+ binding reaction [10, 15]. However, the resultant Hill coefficient was <4 and still lower than the experimental value of 4–7 in the overall Ca2+–F b relationship in the experiments . Thus, we readjusted the number of Ca2+ (nCa) bound to the TS only for a phenomenological description. We found that an assumption of simultaneous binding of five Ca2+ (nCa = 5) gave an overall Hill coefficient of 6–7 in the experimental Ca2+–F b relation (Fig. 2). Thus, the binding rate of Ca2+ to the TS is calculated as,
The concentrations of all states in the scheme of state transition (Fig. 1) are given in terms of TS. Thus, it was assumed that a single TS controls simultaneously a number nCa of CBs. The concentration of the total myosin S1 segment ([mS1]) is given as
The attachment (f) and detachment (g) rates of CB in the NL model are described by Eqs. 4 and 5 , respectively, as a function of the half sarcomere length (halfSL) or the eCB elongation (h w) to reconstruct the transient kinetics evoked by applying rapid changes in length or external mechanical load.
F h in Eq. 5 was defined as described by Eq. 6 in Negroni et al. . The relaxation steps Q 12 and Q 13 (for the Ca2+-unbound conformations) mainly occur during withdrawal of cytosolic Ca2+ and were described using the same rate constant as applied to corresponding state transitions (Q 7 and Q 8) of the Ca2+-bound states, respectively, assuming a competitive binding of ATP in place of ADP.
The rate functions and equations for calculating state transitions are described in the ESM Tables S2 and S3. Equations and parameters to determine F b are the same as those in the original NL model, with the exception of a slight modification of B w, as listed in ESM Table S1.
Concentration of the myosin S1 segment and magnitude of contraction
In the cell models published to date [26,27,28,29,30,31,32,33,34], variable troponin concentrations of approximately 70 μM were assumed. For the [mS1], a molar ratio of actin:troponin concentration of 7:1  and a ratio of myosin S1 segment:actin of 1:4.1  were reported in cardiac muscle. These findings give a ratio of myosin S1 segment:troponin of 1.71, or 119.7 µM [mS1] in the previous cell models. This [mS1] is well within the range of direct biochemical measurements of the myosin head in ventricular tissue in various species, which are 200 (in pig ), 157 (in guinea pig ), and 144 (in rabbit ) µmol/kg wet weight. Here, we tentatively adopt the [mS1] of 200 µM.
In the NL model  the converting factors (stiffness of elastic CB structure) of A p and A w were used to calculate the magnitude of developed tension for a unitary cross-section area (in mm2) of cardiac muscle for the powered and weak CB concentrations (in µM), respectively. In the present study, CB force (F b) is given by the sum of weak CB force (F bw) and powered CB force (F bp), as the following equation (in mN/mm2).
Here, the h p and h w represent elongation of the elastic component of the strong-bound CB states (cAM Ds and AM Ds) and weak-bound CB state (cAM DPw), respectively. Note that variable CB elongations within the whole population of CB within a cell are represented by a single eCB having an ‘average’ CB elongation, h p and h w, for the powered and weak eCBs, respectively in the NL model. The rate of change in eCB elongation and the amount of Ca2+ bound to the TS were calculated in the same way as in the original NL model. The force of the parallel elastic element (F p) is given by Eq. 14.
where hSL 0 is the slack length of the parallel elastic element. The rate of change in the eCB elongation (h w and h p) was calculated essentially in the same way as in the original NL model. The parameter values were adopted from Negroni et al. , as described in ESM Table S1.
The ATP consumption rate is calculated from two state transitions from cAM Ds to cAM T (Q 7) and from AM Ds to AM T (Q 12).
Time-integration of ordinary differential equations
The numerical integration was performed by Euler’s method with a time step of 0.01 or 0.005 ms. The Ca2+ transient was generated by the empirical equations given in the NL model . The model parameters and equations are described in the ESM Tables S1 to S5.
The following dimensions were applied; micromole (μM) for concentrations, milliseconds (ms) for time, milliNewton per millimeter squared (mN/mm2) for force of contraction, and micrometer (µm) for length. All codes of the simulation program were prepared using Microsoft Visual Basic on Visual Studio Community 2013 (Microsoft Corp., Redmond, WA).
ATP hydrolysis is activated by [Ca2+] in the new hybrid contraction model
The AM ATPase is activated indirectly by Ca2+, and repetitive cycles of hydrolysis are maintained in the presence of [Ca2+]. The rate of ATP hydrolysis with accompanying developed tension F b was examined by applying various [Ca2+] to the hybrid model. As the magnitude of F b is dependent on halfSL, the relationship between [Ca2+] and the rate of ATP hydrolysis was examined at different halfSL (Fig. 2a, b). Although the underlying state transitions following Ca2+ binding, including the A 0–A 3 steps are quite complicated (Fig. 1), it would appear that these [Ca2+]–F b or [Ca2+]–v ATPase relationships conform well to a saturation kinetic mechanism (Eq. 15).
To compare with the experimental [Ca2+]–F b relationship obtained by Kentish et al.  and by Gwathmey and Hajjar , we measured the Hill coefficient (n H) by modifying Eq. 15 to a linear form of relationship for F b or ATPase activity (v ATPase) with V max, Fb or V max, ATPase, respectively (Eq. 16; see ESM for the derivation of equations).
Figure 2c and d indeed show approximately linear relationships at every [Ca2+] with nearly constant value of n H. The slope of the relationship was determined as a mean n H over an interval of −0.3 to 0.3 on the abscissa, which increased as 5.348, 5.787, 6.272, and 6.494 at halfSL of 0.8, 0.9, 1.0, and 1.1 µm, respectively. Simultaneously, the K 0.5 decreased as 1.271, 1.105, 0.980, and 0.923 µM in both the [Ca2+]–F b and [Ca2+]–v ATPase relationships. The V max,ATPase increased as 163.03, 245.87, 279.16, and 289.66 µM/ms and the V max,Fb as 737.3, 1112.0, 1262.6, and 1310.1 mN/mm2, at 0.8, 0.9, 1.0, and 1.1 µm halfSL, respectively. Comparison of the relative amplitudes between these values revealed that V max,ATPase is quite proportional to V max,Fb. The slope of the relationship became slightly shallower with increasing [Ca2+], mostly due to the limited number of TS. Note that cooperativity was assumed only for Ca2+ binding to troponin. These results are all quite consistent with the experimental findings [37, 41].
Contraction and ATPase activities evoked by the idealized Ca2+ transient
The proportionality between F b and v ATPase in Fig. 2 is expected in the present hybrid model. Namely, both F b and v ATPase are largely determined by [cAM Ds] and [AM Ds] in both Eqs. 13 and 17. The first term in Eq. 12, which represents the contribution of weakly bound CB, is a minor component of the total F b during the steady presence of Ca2+ because the h w quickly relaxes to a negligibly small h wr (= 0.0001 µm). It may be noted that the denominator in Eq. 17 is constant since [ATP] and [ADP] are both given constants during the time course in Fig. 2.
ATPase activity was examined during the usually developed tension evoked by the transient Ca2+ at a halfSL = 1.05 µm (Fig. 3). The model was activated by the standard transient Ca2+ given by ESM Eqs. S14 and S15. The component of AM T was not included in the calculation of F b, since this state in the new hybrid model represents a sum of (AM T + M T) in the Månsson A 2 and A 3 states. In Fig. 3c, the time course of v ATPase is nearly similar to that of F b in Fig. 3b. During a single stimulus interval of 800 ms, the amount of ATP used was 3969 µM, which gives an average of ATPase turnover rate of 19.84/twitch (=3969/[mS1], [mS1] = 200 µM). The ATP consumption rate via transition step Q 12 during the removal of Ca2+ was much smaller if compared with that in the step Q 7.
Compared to the isometric contraction, the isotonic contraction was characterized by the delayed peak in the developed tension F b and the ATP flux. This delayed peak was caused by the shortening of h p and thereby also of halfSL (Fig. 4e, f), in contrast to the isometric contraction. This decrease in h p largely inverted h w at the onset of contraction, followed by gradual relaxation toward the stable length of h w (0.0001 µm) during contraction. Thus, the CB detachment rate g, given as a function of h w deviation from h wr (Eq. 5), was largely increased (Fig. 4g). On the other hand, the attachment rate f was nearly constant, though increased only slightly. Through these changes in CB kinetics, the development of cAM Ds was largely delayed and its magnitude was smaller, when compared with the isometric contraction. Thereby, the peak tension and the peak of v ATPase were delayed according to the time course of the state cAM Ds. It should be noted that AM Ds becomes a significant component only during the late relaxation phase (Fig. 4).
Taken together, Figs. 3 and 4 indicate that the time course of F b is quite similar to that in the original NL model . Thus, it is evident that the addition of the ATP hydrolysis cycle to the NL model did not affect the general time course of contraction obtained in the new hybrid model.
The length clamp experiment in the hybrid model
The measurements of relaxation time course in F b evoked by applying a step change in halfSL during halfSL-controlled conditions provide essential parameters for the kinetics of the state transitions of the eCB conformations. Thus, we tested our new model by applying the length step shown in Fig. 5a and examined the F b (Fig. 5b) and ATP hydrolysis (Fig. 5c) at a constant [Ca2+] of 0.647 µM. At the onset of the length step, h p decreased to approximately 0 and then rapidly recovered within the subsequent 5 ms due to the rapid intrinsic rate of dh p /dt (Fig. 5e). The negative deflection of h w at the pulse onset caused rapid detachment of the eCB by approximately 80% through the transient increase in g according to Eq. 5. During the following slow relaxation phase of approximately 900-ms interval, the F b recovered due to the re-equilibration of state transitions (Fig. 5d) as in the original NL model.
The time course of v ATPase change (Fig. 5c) was parallel to the time course of [cAM Ds] and [cAM DPw], as shown in Fig. 5d. This finding indicates that ATP consumption should be transiently decreased and at the onset of the shortening length pulses. Thereafter, ATPase activity will gradually recover during the F 4 phase in actual experiments.
The rapid-release experiment in the isotonic contraction mode
The rapid-release experiment revealed the three phases in the time course of halfSL shortening, namely, the initial jump (phase 1), the rapid hyperbolic shortening (phase 2) within the initial 50 ms, and the subsequent slow almost linear shortening (phase 3), as observed typically in skeletal muscle (Fig. 6b, e). These simulation results are quite comparable to the relaxation time course demonstrated in the original NL model . Similar time courses have been obtained in experiments using cardiac muscle preparation, but obviously at a low time resolution due to the intrinsic complexity of the trabeculae in the cardiac tissue .
The force clamp experiment is quite straightforward, namely the product of h p multiplied by the number of eCB is proportional to the applied load (F ext). The quick release instantaneously decreased h p (upper line, Fig. 6d) or negatively deflected h w (lower line, Fig. 6d). The deviation in h w from h wr decreased the number of attached eCBs by accelerating the rate of detachment (g) of eCBs represented by Eq. 5, as in the length jump simulation. This decrease in the number of eCB increased h p to balance the applied F ext. Through these two opposite influences of decreasing and increasing h p, the magnitude of h p relaxed approximately to a new steady level during phases 1–2. The fraction of cAM Ds reached a new equilibrium level through an almost simple sigmoidal time course. Therefore, the rate of ATP hydrolysis was almost linearly related to the level of test external load during the phase 3, as shown in Fig. 6g, which differs from the length clamp experiment.
It should be noted that the magnitude of h p remained nearly equal to 6 nm during the phase 4 shortening. This means that the shortening of halfSL is attributable exclusively to the movement of eCB attachment point along the actin fiber. The rate of shortening in phase 4 was nearly an exponential function of the mechanical load in Fig. 6f, which is in rough agreement with experimental findings (Fig. 6e). On the other hand, the rate of ATP usage increased in proportion to the magnitude of F ext from almost zero in the absence of F ext (Fig. 6g). In summary, the turnover rate of the AM ATPase largely varied depending on both the magnitude of F ext and [Ca2+] in the cytosol ([Ca2+]cyt).
Integration of the new contraction model into the ventricular whole cell model
Parameters of the new hybrid model of contraction were also tuned by integrating the model in the human ventricular cell model (HuVEC model ). This was done because whole cell ATP consumption by the S1 segment has been well established by macroscopic measurements rather than by the study of dissected preparations in vitro. It should be noted that the estimation of ATP consumption by ionic pump activities in the whole cell cardiac cell models is now well established by the refined magnitude of ionic fluxes underlying the membrane excitation [42, 43] as well as by the Ca2+ fluxes underlying the cytosolic Ca2+ dynamics regulated by the sarcoplasmic reticulum (SR) [2, 3, 44].
Figure 7 shows the action potential, the major ionic currents, and the developed tension evoked by the cytosolic Ca2+ transient at a stimulus interval of 800 ms. In the preliminary simulation, we found that the size of the transient Ca2+ of approximately 0.5 μM assumed in most of the human ventricular models [42, 43], including the HuVEC model, failed to evoke the appropriate magnitudes of F b. We therefore readjusted the size of the transient Ca2+ to approximately 0.9 μM by increasing the total SR volume from the 6–8% of total cell volume formerly used and the ratio of the SR releasing site to a total SR from 30 to 40%. These changes did not modify the time course of Ca2+ content in the SR. Under these conditions, ATP consumption by the contraction was 58.6% of the sum of ATP consumption by the myofilaments, SR Ca2+-ATPase (SERCA), Na/K, and plasma membrane Ca2+ ATPase (PMCA) per cycle. This value is quite consistent with the approximately 60% consumption by the AM–ATPase in the beating blood perfused heart . Other components, including SERCA, Na/K-ATPase, and PMCA, used 36, 4.5, and 0.6% of the total ATP consumption, respectively, under the physiological condition. A slightly different result was reported by Schramm et al. , who attributed 76% of the whole ATP consumption to AM–ATPase and, 15 and 9% to the SERCA and Na/K-ATPase, respectively.
Applicability of the present hybrid model to estimate the myocardial ATP usage in the simple blood circulation model
The NL model has been used to examine the hemodynamics in the multi-scale model of the cardiovascular system using ‘Laplace’s left ventricular hemispherical model’ [19, 30, 46]. We replaced the hybrid model for the original NL model in the hemodynamic model of Utaki et al.  and estimated ATP consumption during the pressure–volume (PV) trajectory. As expected from the faithful reproduction of the developed tension described above (Figs. 2, 3), the PV trajectory obtained by the new system model seems to be quite consistent with the results of previous studies (Fig. 8). The pressure–volume area (PVA) of the left ventricle, which correlates well with myocardial oxygen consumption per beat , changed (enlarged) with increasing preload scale (reviewed by Utaki et al. ) (standard K rp × 0.8, × 1.0, × 1.2, × 1.4, × 1.6, and × 1.8 in Fig. 8). ATP usage in the hybrid contraction model was simultaneously plotted in the same figure, but at different scale (Fig. 8a, shown in gray). During the ejection period, the ATP hydrolysis rate decreased because of the decrease in the number of powered CBs. This is different from the gradual increase in ventricular pressure due to the decrease in the radius of the hemisphere Laplace heart. The relationship between PVA (mmHg × ml) and the corresponding integration of ATP usage (μM/ms) was examined in Fig. 8b. The linear relationship between PVA and v ATPase definitely indicates the relevance of using the hybrid model in calculating ATP consumption.
A new cardiac contraction model was developed by combining the biochemical model elaborated by Månsson  with the biophysical NL model. The new model is the successor to the dynamic CB models developed by Lombardi and Piazzesi , Piazzesi et al. , and Edman et al. , all of which are based on tetanus contraction in skeletal muscle. In contrast, the newly proposed NL model of the TS reconstructs F b in the cardiac myocytes over a wide variety of physiological experimental findings. In the new hybrid model, the scheme of state transition among the A 0–A 3 states in the Månsson model was used as it is, but the rate constants for individual state transitions were simplified. In the original Månsson model, the rate was calculated by assuming a Gibbs free energy profile as a function of the distance (x) between the CB head and the nearest binding site on actin fiber. In the hybrid model, a constant rate was used for individual state transitions of eCB by referring to the rate change based on the Gibbs free energy profile. The new model inherited well the major characteristics of both types of the two models, such as the concentration–response relation of [Ca2+]–F b or [Ca2+]–v ATPase (Fig. 2), the time course of the developed tension activated by the intracellular transient Ca2+ (Figs. 3, 4), the turnover rate of the AM ATPase, the F b responses to jumps in the halfSL (Fig. 5), and F ext (Fig. 6) under the isometric and isotonic contraction modes, respectively, and the load–velocity relation (Fig. 6). When the new hybrid model was incorporated into a ventricular cell model, ATP consumption by contraction was approximately 60% of the whole cell ATP usage at a cycle length of 0.8 s (Fig. 7). When the new model was implemented in the hemodynamic blood circulation model, ATP consumption was proportional to the PVA of the hemisphere Laplace heart (Fig. 8).
Simultaneous reconstruction of ATP consumption and the developed tension
The minimum requirement of any mathematical model of cardiac contracting fibers should be the capture of the following four essential steps of the CB cycle tightly coupled to the accompanying ATP hydrolysis:
The binding of ATP in place of ADP to the catalytic domain of the S1 segment rapidly detaches the myosin head from the actin binding site, resulting in relaxation of the rigor state [49,50,51,52].
ATP is hydrolyzed in association with a structural change of a swing of the myosin lever arm (a recovery stroke), leaving products of ATP hydrolysis, ADP, and Pi bound to the active site of myosin head [50, 53, 54].
Ca2+ binding to troponin dislocates the tropomyosin complex during the time transient Ca2+ is within the cell  and thereby allows the myosin head (carrying ADP and Pi) to bind with actin filament, forming a weakly bound state of CB [51, 56].
The attachment of the myosin head to the actin binding site causes an approximately 100-fold increase of the rate of Pi release . Dissociation of Pi from the S1 segment is tightly coupled with the power stroke of the CB, resulting in the sliding motion between the thin and thick filaments to stretch the elastic elements in the CB [57,58,59], which is responsible for the developed tension F b.
Steps 1, 2, and 4 are precisely described in the biochemical models (Månsson model), while the NL model simulated well step 3 based on biophysical experimental findings, such as the time course of F b evoked by the transient Ca2+, the response of the shortening of the halfSL evoked by the rapid release protocol in the isotonic contraction, and the time course of F b recovery evoked by the shortening step of the halfSL in the isometric contraction mode. The hyperbolic increase in shortening velocity with decreasing external load was reconstructed based on the time course of halfSL shortening reconstructed by the numerical integration in the model, by applying the quick release protocol according to the experimental findings. In the Månsson model, this relationship was theoretically predicted indirectly from the transition states. The present model clearly proposes that ATP usage is minimum in the absence of external load and increases in proportion to the magnitude of F ext over the physiological range (Fig. 6g).
We have shown that the new hybrid model successfully calculated the rate of ATP hydrolysis simultaneously with the accompanying development of F b. Thus, the presented hybrid model largely facilitated development of new whole cell models for analyzing cardiac energetics, development of the force of contraction, as well as the ATP dependency of the muscle contraction, as shown in Fig. 7. The merit of using the NL model may be largely due to the introduction of the eCB, which represents the average behavior of the whole population of CBs within a cell. This kind of spatial average model of the whole population of CBs has been used in a variety of simplified models (for references see Introduction). In the NL model, the behavior of the eCB was thoroughly fitted to the classical experimental findings ([60,61,62] and models [8, 63, 64]). The time-dependent change of eCB elongation, h p or h w, is continuous as described in ESM Eqs. S3 or S4, and the rate g of detachment of CB is given as a function of h w (2008 NL model ) or h p (1996 NL model ). This g provides the basis for explaining the halfSL dependence of the [Ca2+]–F b relationship  (Fig. 2). The relationship between the F ext and dh p /dt is also explained by the eCB kinetics, the rate of quick recovery of F b at the onset of the halfSL jump, and the time course of the shortening of halfSL evoked by the release of F ext.
In the original model of Lombardi and Piazzesi , the fractions of CBs in the A l, A 2, and A3 states were calculated by applying the numerical integration method with a discrete interval of (Δx) of <0.5 nm at each x-point. It was also assumed that the elongation of CB was increased by an interval of the neighboring actin binding sites for each state transition of A 1–A 2, and A 2–A 3. When calculating the force–velocity relationship, the smooth continuous change in the velocity as well as the force of CBs were obtained theoretically by averaging for whole population of CB states without reconstructing the experimental response to quick release. Using a linear stiffness of the CB elongation, the force of contraction (for example, in units of mN/mm2) was calculated from the sum of A 2 (weakly bound) + A 3 (strongly bound) over the range of x.
Negroni and Lascano elegantly simplified the calculation of developed tension and the CB kinetics by assuming the eCB might roughly represent the average CB elongation h p or h w . Consequently, the unitary developed tension (f i ) is calculated with a stiffness (a) of the elastic structure (f i = a × h), and the movement rate of the CB head along the actin fiber is calculated as dh/dt. The dh/dt is defined for both h p and h w as:
It should be noted that Eq. 18 is equivalent to ESM Eq. S3 provided that halfSL is constant during the time interval under consideration, as in the L 4 phase in the isotonic contraction (Fig. 6). Thus, dh p /dt is used to calculate the movement of the CB head when attached. In the detached CB head, h relaxes quickly to the resting elongation (h r). Note that the x 0 giving the energy minima in the detached A 3 state of the Månsson model  is about 7–8 nm distant from the corresponding position of energy minima in the A 1 and A 2 state, which is close to h r (=6 nm) in the NL model, although Piazzesi et al.  and Edman et al.  assumed slightly different x 0 of energy minima. In the Månsson model the binding rate is described as a function of the velocity of the CB head movement along the actin fiber ; this assumption may correspond to the h-dependent detaching rate of g in the NL model.
Turnover rate of ATP hydrolysis in the new hybrid model
The sequential steps of Q 7–Q 9 and Q 12–Q 13 in the new hybrid model were adopted from the Månsson model, in which the weakly bound CB was separated from the powered CB according to experimental findings [65,66,67,68,69,70] and then combined with the Q 1–Q 4 and Q 10–Q 11 steps of the 2008 version of the NL model. This modification did not largely modify the time course of developed tension when compared to the NL model because the newly implemented cycle of ATP hydrolysis is much faster than that in the entire state transition cycle in the original NL model. In the new hybrid model, ATP binding to the catalytic domain is also assumed in the unbound Ca2+ state, AM Ds, which appears during the relaxing phase (shown in Fig. 1). This assumption is justified by the recent finding that ATP hydrolysis occurs in the reconstructed myosin fiber in the absence of both Ca2+ and actin . In our simulation (Fig. 3), the contribution of AM Ds, however, was only 1.5% of the total ATP consumption.
The average ATP turnover rate of the new contraction model was 19.85/twitch in Fig. 3. Namely, quite consistent rates of 2.7–3.3/s [71,72,73] have been reported in a rat model of chemically skinned trabeculae in the presence of a saturating concentration of Ca2+ at room temperature (20–21 °C). If these rates are corrected for the physiological temperature using the Q 10 of approximately 3.5 obtained by de Tombe and Stienen , an experimental turnover rate of 20–24/s is expected, which is slightly larger than the present simulation results. However, it should be noted that the rate of ATP consumption depends on several factors, such as the difference between isotonic and the isometric contraction modes (Figs. 3, 4), the [Ca2+]cyt (Fig. 2), the external load (Fig. 5), and the halfSL (Figs. 5, 6). It should also be noted that the measurements of the turnover rate as well as Q 10 showed variations. A much larger ATP turnover rate (approx. 10/s per myosin head at 24 °C) has been reported in skinned rat trabeculae  and, in contrast, a lower value has been described in a pig cardiac preparation (approx. 0.5/s at 24 °C ) and in a rat preparation (approx. 3/s at 20 °C) using different experimental protocols [41, 71,72,73]. The temperature coefficient Q 10 of ATP hydrolysis is also variable (Siemankowski et al. , Q 10 = approx. 5; Burchfield and Rall , Q 10 = approx. 4; de Tombe and Stienen , Q 10 = approx. 3.5). Species-specific differences in the ATP hydrolysis rate might be attributed to differences in the composition of isoforms of the S1 myosin isoform .
The variability of the rate of ATP consumption, as described above, indicates the difficulty in comparing the turnover rate of ATPase activity in simulation and experiments due to the high variability of the recording conditions and data possibly not being completely described in experiments. Thus, any comparison of the relative weight of ATP consumption within a given cell model between the contraction and the ion pumps may be relevant in testing AM ATPase activity. Although additional fine-tuning of several parameters was required in both the contraction model and the cell model, the whole cell model present here reconstructed well the experimental measurement of approximately 60% consumption by the contraction with reference to the sum of contraction and ion pumps (Fig. 7).
The stiffness of single CB in ventricular whole cell
In the present study, the CB concentration of ([CB] = 200 µM) was assumed by referring to corresponding values in conventional cardiac cell models and also to the direct biochemical measurements of CB (see Materials and methods). In the study reported here, we have examined if the [CB] assumed in our model is consistent with the stiffness of single CB of ε = 0.5–3 pN/nm which has been assumed to date in the Månsson model for calculating the F b.
The value of ε is calculated by the following equation using the force of contraction (f b) and the elongation of the single CB (h).
A representative value set was F b = 20 mN/mm2, the fraction of activated CBs f a = 0.05, h = 6 nm, and halfSL = 1 µm. n is the number of CBs within a volume of muscle fiber given as a product of [halfSL × unit cross-section area (A) = 1 mm2] and is determined as:
where, N A is the Avogadro’s number. The ε thus obtained was 0.55 pN/nm, which is well within the range of 0.5–3 pN/nm assumed by Månsson  and Månsson et al. .
Limitation of the new hybrid model
Taken together, the simulation results of the hybrid model presented here are in good agreement to experimental data published to date. However, the limitation of this new contraction mode is obvious, since the ability to generate realistic response does not prove that the underlying biophysical mechanisms are correctly represented.
Suga H (1990) Ventricular energetics. Physiol Rev 70:247–277
Asakura K, Cha CY, Yamaoka H, Horikawa Y, Memida H, Powell T, Amano A, Noma A (2014) EAD and DAD mechanisms analyzed by developing a new human ventricular cell model. Prog Biophys Mol Biol 116:11–24
Himeno Y, Asakura K, Cha CY, Memida H, Powell T, Amano A, Noma A (2015) A human ventricular myocyte model with a refined representation of excitation-contraction coupling. Biophys J 109:415–427
Lombardi V, Piazzesi G (1990) The contractile response during steady lengthening of stimulated frog muscle fibres. J Physiol 431:141–171
Månsson A (2010) Actomyosin-ADP states, interhead cooperativity, and the force–velocity relation of skeletal muscle. Biophys J 98:1237–1246
Månsson A, Rassier D, Tsiavaliaris G (2015) Poorly understood aspects of striated muscle contraction. Biomed Res Int. doi:10.1155/2015/245154
Piazzesi G, Lombardi V (1995) A cross-bridge model that is able to explain mechanical and energetic properties of shortening muscle. Biophys J 68:1966–1979
Landesberg A, Sideman S (1994) Mechanical regulation of cardiac muscle by coupling calcium kinetics with cross-bridge cycling: a dynamic model. Am J Physiol Heart Circ Physiol 267:H779–H795
Razumova MV, Bukatina AE, Campbell KB (1999) Stiffness-distortion sarcomere model for muscle simulation. J Appl Physiol 87:1861–1876
Razumova MV, Bukatina AE, Campbell KB (2000) Different myofilament nearest-neighbor interactions have distinctive effects on contractile behavior. Biophys J 78:3120–3137
Campbell KB, Razumova MV, Kirkpatrick RD, Slinker BK (2001) Nonlinear myofilament regulatory processes affect frequency-dependent muscle fiber stiffness. Biophys J 81:2278–2296
Rice JJ, Wang F, Bers DM, de Tombe PP (2008) Approximate model of cooperative activation and crossbridge cycling in cardiac muscle using ordinary differential equations. Biophys J 95:2368–2390
Negroni JA, Lascano EC (1996) A cardiac muscle model relating sarcomere dynamics to calcium kinetics. J Mol Cell Cardiol 28:915–929
Negroni JA, Lascano EC (2008) Simulation of steady state and transient cardiac muscle response experiments with a Huxley-based contraction model. J Mol Cell Cardiol 45:300–312
Rice JJ, Stolovitzky G, Tu Y, de Tombe PP (2003) Ising model of cardiac thin filament activation with nearest-neighbor cooperative interactions. Biophys J 84:897–909
Hussan J, de Tombe PP, Rice JJ (2006) A spatially detailed myofilament model as a basis for large-scale biological simulations. IBM J Res Dev 50:582–600
Kentish JC (1986) The effects of inorganic phosphate and creatine phosphate on force production in skinned muscles from rat ventricle. J Physiol 370:585–604
de Tombe PP, ter Keurs HE (1991) Sarcomere dynamics in cat cardiac trabeculae. Circ Res 68:588–596
Moriarty TF (1980) The law of Laplace. Its limitation as a relation for diastolic pressure, volume, or wall stress of the left ventricle. Circ Res 46:321–331
Gibbs CL (1978) Cardiac energetics. Physiol Rev 58:174–254
Lab MJ, Allen DG, Orchard CH (1984) The effects of shortening on myoplasmic calcium concentration and on the action potential in mammalian ventricular muscle. Circ Res 55:825–829
Yue DT, Marban E, Wier WG (1986) Relationship between force and intracellular [Ca2+] in tetanized mammalian heart muscle. J Gen Physiol 87:223–242
Janssen PM, Stull LB, Marban E (2002) Myofilament properties comprise the rate-limiting step for cardiac relaxation at body temperature in the rat. Am J Physiol Heart Circ Physiol 282:H499–H507
Stehle R, Iorga B (2010) Kinetics of cardiac sarcomere processes and rate-limiting steps in contraction and relaxation. J Mol Cell Cardiol 48:843–850
Negroni JA, Morotti S, Lascano EC, Gomes AV, Grandi E, Puglisi JL, Bers DM (2015) β-adrenergic effects on cardiac myofilaments and contraction in an integrated rabbit ventricular myocyte model. J Mol Cell Cardiol 81:162–175
Hunter PJ, McCulloch AD, Ter Keurs HEDJ (1998) Modelling the mechanical properties of cardiac muscle. Prog Biophys Mol Biol 69:289–331
Jafri MS, Rice JJ, Winslow RL (1998) Cardiac Ca2+ dynamics: the roles of ryanodine receptor adaptation and sarcoplasmic reticulum load. Biophys J 74:1149–1168
Matsuoka S, Sarai N, Kuratomi S, Ono K, Noma A (2003) Role of individual ionic current systems in ventricular cells hypothesized by a model study. Jpn J Physiol 53:105–123
Matsuoka S, Sarai N, Jo H, Noma A (2004) Simulation of ATP metabolism in cardiac excitation–contraction coupling. Prog Biophys Mol Biol 85:279–299
Coutu P, Metzger JM (2005) Genetic manipulation of calcium-handling proteins in cardiac myocytes. II. Mathematical modeling studies. Am J Physiol Heart Circ Physiol 288:H613–H631
Okada JI, Sugiura S, Nishimura S, Hisada T (2005) Three-dimensional simulation of calcium waves and contraction in cardiomyocytes using the finite element method. Am J Physiol Cell Physiol 288:C510–C522
Shim EB, Amano A, Takahata T, Shimayoshi T, Noma A (2007) The cross-bridge dynamics during ventricular contraction predicted by coupling the cardiac cell model with a circulation model. J Physiol Sci 57:275–285
Shim EB, Jun HM, Leem CH, Matusuoka S, Noma A (2008) A new integrated method for analyzing heart mechanics using a cell-hemodynamics-autonomic nerve control coupled model of the cardiovascular system. Prog Biophys Mol Biol 96:44–59
Soltis AR, Saucerman JJ (2010) Synergy between CaMKII substrates and β-adrenergic signaling in regulation of cardiac myocyte Ca2+ handling. Biophys J 99:2038–2047
Potter JD (1974) The content of troponin, tropomyosin, actin, and myosin in rabbit skeletal muscle myofibrils. Arch Biochem Biophys 162:436–441
Murakami U, Uchida K (1985) Contents of myofibrillar proteins in cardiac, skeletal, and smooth muscles. J Biochem 98:187–197
Kuhn HJ, Bletz C, Rüegg JC (1990) Stretch-induced increase in the Ca2+ sensitivity of myofibrillar ATPase activity in skinned fibres from pig ventricles. Pflug Arch 415:741–746
Barsotti RJ, Ferenczi MA (1988) Kinetics of ATP hydrolysis and tension production in skinned cardiac muscle of the guinea pig. J Biol Chem 263:16750–16756
Hooijman P, Stewart MA, Cooke R (2011) A new state of cardiac myosin with very slow ATP turnover: a potential cardioprotective mechanism in the heart. Biophys J 100:1969–1976
Gwathmey JK, Hajjar RJ (1990) Relation between steady-state force and intracellular [Ca2+] in intact human myocardium. Index of myofibrillar responsiveness to Ca2+. Circulation 82:1266–1278
de Tombe PP, Stienen GJM (1995) Protein kinase A does not alter economy of force maintenance in skinned rat cardiac trabeculae. Circ Res 76:734–741
Grandi E, Pasqualini FS, Bers DM (2010) A novel computational model of the human ventricular action potential and Ca transient. J Mol Cell Cardiol 48:112–121
O’Hara T, Virág L, Varró A, Rudy Y (2011) Simulation of the undiseased human cardiac ventricular action potential: model formulation and experimental validation. PLoS Comput Biol 7:e1002061
Hinch R, Greenstein JL, Tanskanen AJ, Xu L, Winslow RL (2004) A simplified local control model of calcium-induced calcium release in cardiac ventricular myocytes. Biophys J 87:3723–3736
Schramm M, Klieber HG, Daut J (1994) The energy expenditure of actomyosin-ATPase, Ca2+-ATPase and Na+, K+-ATPase in guinea-pig cardiac ventricular muscle. J Physiol 481:647–662
Utaki H, Taniguchi K, Konishi H, Himeno Y, Amano A (2016) A method for determining scale parameters in a hemodynamic model incorporating cardiac cellular contraction model. Adv Biomed Eng 5:32–42
Piazzesi G, Francini F, Linari M, Lombardi V (1992) Tension transients during steady lengthening of tetanized muscle fibres of the frog. J Physiol 445:659–711
Edman KAP, Månsson A, Caputo C (1997) The biphasic force–velocity relationship in frog muscle fibres and its evaluation in terms of cross-bridge function. J Physiol 503:141–156
Siemankowski RF, White HD (1984) Kinetics of the interaction between actin, ADP, and cardiac myosin-S1. J Biol Chem 259:5045–5053
Rayment I, Holden HM, Whittaker M, Yohn CB, Lorenz M, Holmes KC, Milligan RA (1993) Structure of the actin–myosin complex and its implications for muscle contraction. Science 261:58–65
Spudich JA (1994) How molecular motors work. Nature 372:515–518
Yount RG, Lawson D, Rayment I (1995) Is myosin a “back door” enzyme? Biophys J 68:44S–49S
Lymn RW, Taylor EW (1971) Mechanism of adenosine triphosphate hydrolysis by actomyosin. Biochemistry 10:4617–4624
Sugi H, Minoda H, Inayoshi Y, Yumoto F, Miyakawa T, Miyauchi Y, Tanokura M, Akimoto T, Kobayashi T, Chaen S, Sugiura S (2008) Direct demonstration of the cross-bridge recovery stroke in muscle thick filaments in aqueous solution by using the hydration chamber. Proc Natl Acad Sci USA 105:17396–17401
Ford LE (1991) Mechanical manifestations of activation in cardiac muscle. Circ Res 68:621–637
Shimizu H, Fujita T, Ishiwata SI (1992) Regulation of tension development by MgADP and Pi without Ca2+. Role in spontaneous tension oscillation of skeletal muscle. Biophys J 61:1087–1098
Trentham DR, Eccleston JF, Bagshaw CR (1976) Kinetic analysis of ATPase mechanisms. Q Rev Biophys 9:217–281
Geeves MA, Goody RS, Gutfreund H (1984) Kinetics of acto-S1 interaction as a guide to a model for the crossbridge cycle. J Muscle Res Cell Motil 5:351–361
Goldman YE, Brenner B (1987) Special topic: molecular mechanism of muscle contraction. Annu Rev Physiol 49:629–636
Huxley AF, Simmons RM (1971) Proposed mechanism of force generation in striated muscle. Nature 233:533–538
Huxley AF (1974) Muscular contraction. J Physiol 243:1–43
Ford LE, Huxley AF, Simmons RM (1977) Tension responses to sudden length change in stimulated frog muscle fibres near slack length. J Physiol 269:441–515
Izakov VY, Katsnelson LB, Blyakhman FA, Markhasin VS, Shklyar TF (1991) Cooperative effects due to calcium binding by troponin and their consequences for contraction and relaxation of cardiac muscle under various conditions of mechanical loading. Circ Res 69:1171–1184
Peterson JN, Hunter WC, Berman MR (1991) Estimated time course of Ca2+ bound to troponin C during relaxation in isolated cardiac muscle. Am J Physiol Heart Circ Physiol 260:H1013–H1024
Brenner B, Schoenberg M, Chalovich JM, Greene LE, Eisenberg E (1982) Evidence for cross-bridge attachment in relaxed muscle at low ionic strength. Proc Natl Acad Sci USA 79:7288–7291
Geeves MA, Holmes KC (1999) Structural mechanism of muscle contraction. Annu Rev Biochem 68:687–728
Baumann BA, Liang H, Sale K, Hambly BD, Fajer PG (2004) Myosin regulatory domain orientation in skeletal muscle fibers: application of novel electron paramagnetic resonance spectral decomposition and molecular modeling methods. Biophys J 86:3030–3041
Ferenczi MA, Bershitsky SY, Koubassova N, Siththanandan V, Helsby WI, Panine P, Roessle M, Narayanan T, Tsaturyan AK (2005) The “roll and lock” mechanism of force generation in muscle. Structure 13:131–141
Wu S, Liu J, Reedy MC, Tregear RT, Winkler H, Franzini-Armstrong C, Sasaki H, Lucaveche C, Goldman YE, Reedy MK, Taylor KA (2010) Electron tomography of cryofixed, isometrically contracting insect flight muscle reveals novel actin–myosin interactions. PLoS One 5:e12643
Smith DA, Geeves MA, Sleep J, Mijailovich SM (2008) Towards a unified theory of muscle contraction. I: foundations. Ann Biomed Eng 36:1624–1640
Kentish JC, Stienen GJ (1994) Differential effects of length on maximum force production and myofibrillar ATPase activity in rat skinned cardiac muscle. J Physiol 475:175–184
Ebus JP, Stienen GJ, Elzinga G (1994) Influence of phosphate and pH on myofibrillar ATPase activity and force in skinned cardiac trabeculae from rat. J Physiol 476:501–516
Ebus JP, Stienen GJM (1996) ATPase activity and force production in skinned rat cardiac muscle under isometric and dynamic conditions. J Mol Cell Cardiol 28:1747–1757
de Tombe PP, Stienen GJM (2007) Impact of temperature on cross-bridge cycling kinetics in rat myocardium. J Physiol 584:591–600
Wannenburg T, Janssen PM, Fan D, de Tombe PP (1997) The Frank-Starling mechanism is not mediated by changes in rate of cross-bridge detachment. Am J Physiol Heart Circ Physiol 273:H2428–H2435
Siemankowski RF, Wiseman MO, White HD (1985) ADP dissociation from actomyosin subfragment 1 is sufficiently slow to limit the unloaded shortening velocity in vertebrate muscle. Proc Natl Acad Sci USA 82:658–662
Burchfield DM, Rall JA (1986) Temperature dependence of the crossbridge cycle during unloaded shortening and maximum isometric tetanus in frog skeletal muscle. J Muscle Res Cell Motil 7:320–326
Stienen GJ, Kiers JL, Bottinelli R, Reggiani C (1996) Myofibrillar ATPase activity in skinned human skeletal muscle fibres: fibre type and temperature dependence. J Physiol 493:299–307
The authors received no external funding for this research.
Conflict of interest
The authors declare that they have no conflict of interest.
This article does not contain any studies with human participants or animals performed by any of the authors.
Electronic supplementary material
Below is the link to the electronic supplementary material.
Rights and permissions
About this article
Cite this article
Muangkram, Y., Noma, A. & Amano, A. A new myofilament contraction model with ATP consumption for ventricular cell model. J Physiol Sci 68, 541–554 (2018). https://doi.org/10.1007/s12576-017-0560-x
- Myofilament model
- Crossbridge kinetics
- Troponin system