Hypothesis: Single Actomyosin Properties Account for Ensemble Behavior in Active Muscle Shortening and Isometric Contraction

Muscle contraction results from cyclic interactions between myosin II motors and actin with two sets of proteins organized in overlapping thick and thin filaments, respectively, in a nearly crystalline lattice in a muscle sarcomere. However, a sarcomere contains a huge number of other proteins, some with important roles in muscle contraction. In particular, these include thin filament proteins, troponin and tropomyosin; thick filament proteins, myosin binding protein C; and the elastic protein, titin, that connects the thin and thick filaments. Furthermore, the order and 3D organization of the myofilament lattice may be important per se for contractile function. It is possible to model muscle contraction based on actin and myosin alone with properties derived in studies using single molecules and biochemical solution kinetics. It is also possible to reproduce several features of muscle contraction in experiments using only isolated actin and myosin, arguing against the importance of order and accessory proteins. Therefore, in this paper, it is hypothesized that “single molecule actomyosin properties account for the contractile properties of a half sarcomere during shortening and isometric contraction at almost saturating Ca concentrations”. In this paper, existing evidence for and against this hypothesis is reviewed and new modeling results to support the arguments are presented. Finally, further experimental tests are proposed, which if they corroborate, at least approximately, the hypothesis, should significantly benefit future effective analysis of a range of experimental studies, as well as drug discovery efforts.


Supporting Text
The parameter values used in the model simulations in the main paper are given in Tables S1 and S2 and the rate functions are given in Table S3. The latter functions are inserted into the ordinary differential equations used to obtain steady-state probability distributions vs the actin-myosin strain coordinate x for the different cross-bridge states 1, 2 . These probability distributions are then used to obtain steady-state force, stiffness, ATP turnover rates etc. 1, 3 . The rate functions in Table  S3 with parameter values as in Tables S1 and S2 are also used to obtain stochastic waiting times between subsequent state transitions and stochastic selection of the actual transition for each event according to the Gillespie algorithm 4 (see further 5 ), when simulating transient events and events involving few (<100) myosin molecules.
The origin of the parameter values are indicated by references in Tables S1 and S2 but are commented on in greater detail below. Importantly, all parameter values are obtained under as similar experimental conditions as possible using isolated myosin (full length, subfragment 1 or heavy meryomyosin) from fast skeletal muscle of the rabbit. Further, all values refer to ionic strengths of 100-200 mM, pH 7-8 and temperature of 30 o C unless otherwise stated in Tables S1-S2.
Parameter values 1: x-values for positions of free energy minima of different states First, a coordinate system that describes the actomyosin distortion (strain) is defined, so that x=0 nm when the free energy of the rigor (AM) state is at its minimum value, thus fixing x3 at x3=0 nm. Moreover, the total power-stroke distance (x11 -x3 = (x11 -x2) + (x2 -x3); see main Fig. 3B), assumed here and in most other studies to involve more than one stroke 6 , has been found to be in the range 7 -13 nm 6 suggesting that x11 is in this range. In the simulations, I use a crossbridge stiffness value (ks=2.8 pN/nm) that I judge to have credible support in single molecule mechanics measurements using full length myosin 7 . This implies a total power-stroke distance (x11-x3) close to the lower bound of the 7-13 nm range because the elastic energy (ks(x1-x3) 2 /2) reaches the free energy of ATP turnover (25 kBT) for x1=8.4 nm. Together with a maximum thermodynamic efficiency of muscle of less than 50 % 8 it follows from this result that (x11-x3) is likely to be lower than 8.4 nm. In our initial modelling 3 we used x11=8.2 nm, later optimized (within experimental uncertainties) in attempts to account for effects of varied concentrations of inorganic phosphate (Pi), to x11=6.7 nm. It should be noted that the latter value relies on the assumption that cross-bridge stiffness is 2.8 pN/nm (see above). If a stiffness value of 1.5 pN/nm is instead assumed, as suggested by some studies 9, 10 , the quantity ks(x11-x3) 2 /2) attains 25 kBT for x11=11.5 nm rather than for x11=8.4 nm, consistent with the rather broad range of the power-stroke distance in the literature (reviewed in 6 ). Now, there is ample evidence for more than one force-generating power-stroke. Whereas some authors have assumed two or more power-strokes of similar amplitudes (e.g. [11][12][13] ), most studies (e.g. [14][15][16] ; reviewed in 6 ) suggest that a large major power-stroke (between the AMDL and AMDH states in present model; amplitude (x11-x2)) is followed by a second smaller stroke, presumably associated with strain-dependent release of ADP (between AMDH and AMD/AM; amplitude x2-x3) (see further main Fig. 3). Single molecule mechanics data suggest that the amplitude x2-x3 is in the range 0.9 -2.5 nm 14,15 . In my initial modelling efforts (when the higher estimate of 2.5 nm 15 was not available) I used x2=1 nm and now continue to do so. However, if a higher numerical value of x2 is corroborated in future studies the model will be changed accordingly. Finally, the parameter value x1 that gives the position of the free energy minimum of the pre-power-stroke state was introduced (in addition to the parameter represented by x11) in efforts to model effects of the small molecular myosin inhibitor blebbistatin 3 . Initially we assumed that x1=x11+1nm but this was later 17 changed to x1=x11+0.5nm to accommodate effects of varied [Pi]. Such changes are within the experimental uncertainties because no real quantitative estimate exists for the difference x1-x11. However, values of about 1 nm or less are reasonable in view of the evidence from X-ray crystallography and arguments that only a small lever arm movement should occur prior to the main powerstroke in order to optimize force generation 18 .
Parameter values 2: differences in free energy minima between cross-bridge states Next in Table S1, I consider differences in minimum free energy levels between neighboring states in the cross-bridge model in main Fig. 3. Measurements in reference 19 at 10 o C and 20 o C, suggest that a quantity corresponding to the sum  Gon+  GPs is in the range of 1.5 -2.8 kBT. The derivation of this range assumes Q10 in the range, 2.7-3.7 19,20 for the attachment rate constant with negligible temperture sensitivity of the reverse rate constant. After some fine-tuning in previous studies, we use the value 1.7 kBT for the sum  Gon+  GPs tentatively assuming that  Gon=0.7 kBT and  GPs =1 kBT based on reasoning in 18 which suggests that both these actomyosin binding strengths are attributed mainly to weak electrostatic interactions. Importantly, the uncertainties in the distribution of the total free energy change (  Gon+  GPs) between  Gon and  GPs are likely to be of negligible significance under physiological conditions. This follows because models that merge the two transitions into one give quite similar predictions as models were the transitions are treated separately 17 . The free energy drop, GLH, associated with the main power-stroke transition (from AMDL to AMDH) is likely to be in the range 11 -19 kBT based on a sum of GLH + GHR in the range 12-20 kBT 6, 7 estimated from power-stroke distances and stiffness measurements from optical tweezers studies of full length myosin. Furthermore, an observed maximum isometric force in optical traps of 17 pN 21 provides a direct estimate of the difference between the free energy minima of the AMDL and AMDH states. With a stiffness of 2.8 pN/nm, a force of 17 pN corresponds to a strain of 6.1 nm and an thus an elastic energy of 2.8 x 6.1 2 /2= 52.1 pN nm ≈ 13 kBT which can be taken as a direct estimate of GLH. Based on these results, together with thermodynamic efficiency of less than 50 % 8 , a quantity x2-x3=1 nm and GHR = 2kBT, we take GLH = 14 kBT in our simulations. The value GHR = 2kBT relies on a value >1.6 kBT for free-energy differences betweeen the two actomyosin-ADP states found in solution studies 22 . Furthermore, a cross-bridge stiffness as 2.8 pN/nm and a quantity x2-x3=2.5 nm (the highest literature value) suggest a free-energy drop GHR in connection with the second power-stroke of 2.1 kBT. In view of these results we use GHR = 2 kBT in the simulations. This value also turns out to be of correct magnitude to account for the high-force deviation of the force-velocity relationship and the molecular effects of the small molecular compound amrinone 16 .
The free energy change associated with the actual Pi-release step (Table S1) is 3.1 kBT, assuming an intracellular Pi-concentration of 0.5 mM. Further, the free energy change associated with the ATP hydrolysis/recovery stroke is kBT ln(10)≈2.3 kBT, based on the equilibrium constant K3=10 (Table S2). Now, adding these values (2.3+3.1) kBT =5.4 kBT to the sum (  Gon+  GPs +GLH + GHR) = 17.7 kBT we arrive at a total sum of 23.1 kBT leaving approximately a 2 kBT drop in free energy to be associated with ATP induced actomyosin dissociation to add up to 25 kBT drop in free energy upon turnover of one ATP molecule. Such a low drop in free energy in the latter step is consistent with a high thermodynamic efficiency (cf. 23 ).

Parameter values 3: rate constants
Starting by considering the rate constants for the ATP hydrolysis (Table S2) on the active site of myosin lumped together with the recovery stroke, litterature data suggest values for the sum of the forward and reverse transitions (k+3+k-3) in the range 300-500 s -1 with an equilibrium constant 24 K3 of 4-10. These values are based on solution kinetics results at lower temperature 24, 25 so they also rely on the assumption of a Q10 value in the range 3-4 from other studies 24 . What is actually used in the simulations is k+3+k-3 = 220 -1 and K3 =10. The low value of k+3+k-3 used, reflects a previously used range of temperatures from 25-30 o C. For the purpose of the present work k+3+k-3 has been kept at its previous value 5,17,[26][27][28] . However, for future work it should be increased about two-fold to the proper range which actually would be expected to further increase the model power-output slightly. Now, assuming k+3=400 s -1 (as appropriate for 30 o C according to above discussion) and that Q10 is 2.8-5 20, 24 for the maximum actomyosin ATP turover rate (Vmax) in solution, then measurements of Vmax at 20-25 o C giving values in the range 13-50 s -1 20, 24 suggest that the parameter kon´, determining the rate constant of cross-bridge attachment, would be in the range 40-150 s -1 . This range would be broader (50-240 s -1 ) if k+3=200 s -1 as actually assumed here. In the simulations, I thus set kon´ = 130 s -1 , approximately in the middle of the latter range.
The constant kP+´, that is a main determinant of the transition rate from the pre-power-stroke state (AMDPP) to the phosphate release state (AMDPiR) (cf. main Fig. 3), was introduced 3 , based on previous structural evidence 18 in modelling to account for the effects of the small molecular myosin inhibitor blebbistatin. The exact numerical value of kP+´ is unknown. However, the previous work 3 suggests that the numerical value is sufficiently high to have limited effect on the ATP turnover rate and on the maximum sliding velocity under physiological conditions (in the absence of drugs). Whereas kP+´ = 1000 s -1 was initially used, an increase of kP+´ from 1000 s -1 to 3000 s -1 was implemented more recently 17 because it dramatically increases the maximum velocity (30 %) with minimal (<3 %) effects on maximum isometric tension, rate of rise of isometric force, steady-state ATP turnover rate and shape (curvature; e.g. maximum power) of the force-velocity relationship. Such an increase in velocity was found essential when other minor changes (in x1 and x11) were implemented to account for effects of varied [Pi] on the force-velocity relationship.
The rate constant of cross-bridge detachment at the end of the power-stroke and in the dragstroke region 3,26,29 is the primary determinant of the maximal shortening velocity although this is also modulated by the value of kP+´ and the degree of linearity of the cross-bridge stiffness as outlined above and in the main paper.
The cross-bridge detachment, with transition from the AMDH state to the MT state in Fig. 3, is governed by a rate constant kdiss: This means that k6, K1 as well as the constants determining the numerical values of k5(x) (primarily k-5) and k2(x) (primarily k2 and xcrit) will have major role in determining the maximum velocity of shortening while also having effects on the shape of the force-velocity relationship. Fortunately, quantiative estimates of k6, k2 and K1 can be found in the litterature 30 for fast rabbit myosin at ionic strengths in the range 100 -200 mM (as used here) and temperature of 25 o C. I directly use the values of k6 and K1 from that paper because they exhibited limited temperature dependence. However, the value of k2 increased with temperature which I take into account here in setting k2=2000 s -1 in the simulations. From the values of k6, k2, and K1 from 30 it follows that k-5 must have a similar magnitude as k2 and k6 to account for the high maximum velocity of shortening. However, at temperatures >25 o C both k2 and rate constants associated with ADP release are likely to be of similar magnitude 30 . Here, I therefore set k-5 = 2000 s -1 , similar to the value of k2 . The parameter xcrit, which defines the strain-dependence of the rate function k2(x) had to be taken as 0.6 nm to account for the observed maximum velocity of shortening in the range 13000-18000 nm per half-sarcomere per second for fast mammalian muscle at 30 o C (see 26 and references therein) with other parameter values as in Tables S1-S2. This was necessary, despite the fact that ultrafast force clamp records using fast myosin subfragment 1 from the mouse at 20 o C suggest that xcrit<0.2 nm. The discrepancy could have different causes. One possibility is that the experimental data from mouse subfragment 1 at 20 o C do not reflect the properties of full length rabbit psoas muscle at 30 o C. Another possibility is that the cross-bridge elasticity is nonlinear also in muscle as found previously for isolated myosin molecules 7 . In the latter case, the high maximum velocity of shortening of mammalian muscle is explained without any strain dependence of k2(x), i.e. with xcrit=0 nm.   kon´ exp(  Gon -ks (x-x1) 2 / kBT + ksw(x-x1) 2 / kBT) Reversal of kon(x) kon´ exp(ks (x-x1) 2 /kBT ksw(x-x1) 2 Table S3 a For further definition of symbols, see main Fig. 3 and for parameter values used, see Tables S1 and S2. More details in 3 .