Mechanosensing in Myosin Filament Solves a 60 Years Old Conflict in Skeletal Muscle Modeling between High Power Output and Slow Rise in Tension

Almost 60 years ago Andrew Huxley with his seminal paper (Huxley, 1957) laid the foundation of modern muscle modeling, linking chemical to mechanical events. He described mechanics and energetics of muscle contraction through the cyclical attachment and detachment of myosin motors to the actin filament with ad-hoc assumptions on the dependence of the rate constants on the strain of the myosin motors. That relatively simple hypothesis is still present in recent models, even though with several modifications to adapt the model to the different experimental constraints which became subsequently available. However, already in that paper, one controversial aspect of the model became clear. Relatively high attachment and detachment rates of myosin to the actin filament were needed to simulate the high power output at intermediate velocity of shortening. However, these rates were incompatible with the relatively slow rise in tension upon activation, despite the rise should be generated by the same rate functions. This discrepancy has not been fully solved till today, despite several hypotheses have been forwarded to reconcile the two aspects. Here, using a conventional muscle model, we show that the recently revealed mechanosensing mechanism of recruitment of myosin motors (Linari et al., 2015) can solve this long standing problem without any further ad-hoc hypotheses.


Introduction
Muscle contraction is generated by the cyclical interaction of myosin motors with the actin filament. Following attachment the myosin motor produces force and relative filaments sliding. This cyclical interaction has been extensively studied and characterized leading to the definition of the myosin-actin cycle. The kinetics of the transitions between the myosin states in the cycle are responsible for both the isometric contraction, when, following activation, the force rises to a maximum steady value, and the isotonic contraction, when the muscle shortens at a constant velocity which increases with the reduction of the load, according to the force-velocity relation Hill (1938). Muscle models based on relatively simple myosin-actin cycles, show that the kinetics needed to simulate the high-power output at intermediate shortening velocities, imply a rate of force development in isometric contraction much faster than that experimentally observed Piazzesi und Lombardi (1995); Huxley und Tideswell (1997); Mansson (2010) ;Caremani u. a. (2013); Marcucci und Yanagida (2012). Despite decades of experimental discoveries and characterization of a more rich environment in the cross bridge cycle, apparently no simple explanation is able to reconcile this controversial behavior. Different additional mechanisms have been introduced on the myosin-actin cycle to solve the problem: either a non-conventional path in cross-bridge cycle Piazzesi und Lombardi (1995), or an active role of the second myosin head Huxley und Tideswell (1997), or a dependence of the rate of attachment on the shortening velocity Mansson (2010), or the possibility that one Caremani u. a. (2013) or more actin monomers become available for the same attached myosin during one ATPase cycle Marcucci und Yanagida (2012). None of these hypotheses, however, provided a definite convincing solution of the problem.
A recent breakthrough on thick filament structure has substantiated the existence of two distinct configurations of the myosin motors in the detached state. In one configuration, motors at rest lie along the thick filament backbone, regularly packed in quasi-helical grooves, as originally demonstrated with cryo-electromicroscopy Woodhead u. a. (2005), while in the second configuration they diffuse towards the thin filament. The presence in relaxed skinned fibers of two populations of myosin motors, one with a low ATP hydrolysis rate (super-relaxed or OFF state) and the other, more disordered (active or ON state) with ten times higher ATPase rate, has been first shown by Cooke in relaxed skinned fibers Stewart u. a. (2010). The muscle activation leads to a transition of the myosin motors from the OFF state to the ON state. Only in this latter state, myosins can bind the actin filaments and perform the power stroke that leads to generation of force and/or filament displacement. Recently published evidence on intact frog muscle fibers at low temperature, shows that the rise of filament stress accompanying high load contraction induces progressive recruitment of motors for the interaction with the actin Linari u. a. (2015). At zero or very low force (less than one tenth of the maximum isometric force T 0 ), only a very few myosin motors are in the ON state, i.e. ready to bind actin filament, and are responsible for the rapid development of the ability to shorten at the maximum velocity Lombardi und Menchetti (1984). The number of motors in the ON state increases during the rise of isometric force and reaches the maximum value at about 0.5 T 0 . In support to the stress dependent transition, keeping to zero the tension of a contracting muscle fiber with the imposition of shortening at the maximum velocity V 0 for 20 − 40 milliseconds, sensibly reduces the population in the active state, favoring the transition from the ON to the OFF state. The change in the proportion between the two populations, affects the time required for the recovery of a steady tension when the isometric condition is restored. The time constant for the rise in tension after initial activation, thus starting from almost all myosins in the OFF state, has been estimated at about τ R = 34 ms. After an unloaded period of 40 ms, corresponding to a shortening of 10% of the muscle length, the tension recovery is faster, with a time constant of τ 10 = 28 ms, because only a part of the ON myosins turns OFF. After an unloaded period of 20 ms, corresponding to a shortening of 5% of the muscle length, the time constant is even lower, τ 5 = 24 ms, because an even minor part of myosins became OFF Linari u. a. (2015). In our view, this imply that the characteristic time scale of the ON-OFF transition is, then, in the same order of magnitude of the attachment detachment process. Consequently, this force regulator may be responsible for the relatively slower rise after activation compared to the relatively faster attachment and detachment rates required for the high power output at intermediate velocities of contraction.
In the present work we test the above hypothesis introducing the experimentally observed mechanosensing mechanism into a recently designed, conventional model of one of the authors Marcucci u. a. (2016). At the best of our knowledge, this is the first time that such a mechanism is embedded in the cross-bridge cycle governing a muscle contraction model. The previous model is able to reproduce a number of experimental data, but it shares the common limit of mathematical models described before: the maximum power output is relatively small compared to the experimental data, because of a relatively high inflection in the force velocity curve ( Fig. 4 in Marcucci u. a. (2016)). Starting from that conventional model, we imposed the observed thick filament tension dependence on the rate constants between the two non-force generating states, here described in term of OFF state and ON state. This modification introduces the mechanosensing mechanism. The description of all other stable states are kept the same. We then compare the simulated behavior to highlight the physiological meaning of the observed property.

Material and Methods
In this paper we test the above mentioned hypothesis comparing the two mathematical models using a Monte-Carlo numerical method. We will refer to the old model as conventional model and to the new one as mechanosensing (MS) model. The conventional model is described in Marcucci u. a. (2016), we resume its main features here and describe the modifications to include the mechanosensing mechanism. The model (see Figure 1) describes a single half-sarcomere with N f il couples of thin and thick filaments and N xb myosin motors per thick filament (compatible with the physiological values in the 2-D simplification, see Table S1 and SI). We simulate the tension-time and force-velocity curves of the whole fiber, under the hypothesis of homogeneous behavior in the series of sarcomeres. Importantly, the experimental reference is given by the frog muscle fiber contractile behavior at 4 • C as described in Piazzesi u. a. (2007) ;Linari u. a. (2015); Lombardi u. a. (1992).
Both filaments are rigid. Actin filaments are constrained to have the same value in the z-direction through the Z-line which represents the coupling term in the mathematical model. z = 0 refers to the optimal length of the sarcomere. Myosin motors are attached through an elastic element to the thick filament backbone and can cyclically interact with an actin filament. Their position is defined by x i (i=1,2...N xb ), the stretch of the motor elastic element. Myosin motor has then five stable states in both models: two non-force-generating states and three force-generating states (including a two step power stroke). In MS model, the non-force-generating states correspond to the super-relaxed, of OFF state, with a negligible ATP hydrolysis, and to the disordered relaxed, or ON, state. In the conventional model, they are representative of the detached and weakly attached state. For simplicity we keep the OFF and ON nomenclature. Motors in the ON state are ready to attach to actin filament forming the strongly attached states. In the conventional model, the OFF and ON transition rates, k 01 and k 10 , respectively, include the cooperativity between myosin motors. Cooperativity affects the transition rates k 01 and k 10 through a parameter γ n , where n is 0 if there are no nearest neighbors and 1 or 2 if one or two, respectively, myosin motors are already weakly or strongly attached. Then k 01 = k 0 01 γ n and k 10 = k 0 10 γ −n Rice u. a. (2003). Such cooperativity is related to the displacement of the tropomyosin filament generated by an attached myosin motor, thus it is related to the thin filament. In MS model, we abandon the cooperativity mechanism, instead we introduce the novelty of the mechanosensing mechanism, i.e. a stress dependent regulation based on the thick filament. The OFF to ON transition rate depends on the force born by the thick filament as observed experimentally Linari u. a. (2015). As described in the introduction, the probability of switching between ON and OFF depends on the tension in the thick filament. Experimentally, the mechanosensing mechanism has been studied through the X-ray diffraction method, in particular analyzing the intensity of M6 reflection at different tensions acting on the thick backbone. It is not straightforward to make a direct relationship between that observable and the myosins distribution in the model since it would require a detailed description of the three-dimensional geometry as well as some hypotheses on the backbone elasticity which may be not linear. Then, as the simplest hypothesis, the OFF to ON rate (k 01 ) is a constant multiplied by a factor ON f (T ), which is a function of the tension borne by the thick filament.  Table S1).
The inverse transition rate, k 10 , has been chosen constant as the simplest Figure 2. ON factor dependence from the relative tension on the thick filament. Experimental data refers to SM6 signal and a direct comparison with the model prediction is not possible (see text). Then, the ON factor has been defined to (i) reach about 5% of ON motors a low forces (less than 0.1 T 0 ), and (ii) saturating at forces higher than 0.5 T 0 .
hypothesis. The ratio of these two rates is chosen to allow for about 5% of myosin motors in the ON state when muscle is relaxed, as hypothesized in the experimental work.
In the ON state, myosins are in equilibrium with the strongly attached state. Both the conventional and MS model have three strongly attached states (S), in either pre-power stroke (S0), first (S1) or second (S2) post-power stroke. In the MS model, myosin motors detach from the strongly attached state switching to the ON state, while in the conventional model the detachment lead to a OFF state. In the OFF state, myosin motors are completely prevented to strongly attach to the actin filament.
From here on the models are equal. The attachment occurs only in the pre-power stroke state, while the detachment can occur in both the pre and the post power stroke state. The attachment to and detachment form the S states follow the classical H57 hypothesis HUXLEY (1957). The attachment rate k 12 is non-zero only in the stretched configuration of the elastic element (x > 0), increases linearly with x up to x lim . The detachment rate, k + 20 or k + 21 for the conventional (S to OFF) or the MS (S to ON) model respectively, also increases linearly with x > 0, and attains a very high value when x < 0 (k − 20 or k − 21 ). The former hypothesis, related to the Brownian search and catch mechanism, has an original observations in Myosin V Fujita u. a. (2012), and has been extensively used in muscle modeling under different modifications. If the absolute value of x exceeds 20 nm, mechanical dislodging occurs. Numerical values used for the models are resumed in Table S1.
In both models Calcium concentration is saturated and TT units reach the full activation within the first millisecond Fusi u. a. (2014) (see SI in Marcucci u. a. (2016)). The strongly attached state is described as a continuous energy landscape with three minima, corresponding to a pre-power stroke and a power stroke with two steps. The central region of the energy landscape is defined as: x, as said, is the stretch of the elastic element in each myosin motor, respect its anchor on the thick filament. x 0 is the stretch at the time of attachment, which depends on the H57 hypothesis plus a random component within [-2.5 nm, 2.5 nm] to simulate the actin diameter Huxley und Tideswell (1996). ATP energy release favors the power stroke by biasing the biochemical energy toward the post power stroke. To simulate this bias, we added a linear drop F atp of 8KT (K is the Boltzmann constant and T denotes the absolute temperature) every d, the distance between two consecutive minima, to the flat sinusoidal part. α d is a constant angle that adjusts for the constant ATP drop by shifting the first minimum to x = x 0 . Meanwhile, the convex part is expressed by a polynomial, ensuring a continuous first derivative. The energy barrier in the sinusoidal function, H, is chosen to be 5.8KT . We can then approximate the minimum in the energy as a parabola with stiffness 48 pN/nm, much higher than the 2 pN/nm stiffness of the myosin motor. The lengths of the half-actin and myosin filaments are LA=1224 nm and LB=925 nm, respectively. The myosin filament includes LB=50 nm of bare zone with no myosin motors. We included 76 myosin acting on one thin filament.

Results and Discussion
Parameters are defined following some assumptions required to reproduce the contractile response of frog muscle fibers at low temperature. Some assumptions are shared between the two models: the time constant of the rising phase after activation must be comparable to the one observed experimentally in Linari u. a. (2015) (τ R = 34 ms); about one third of myosin motors are attached during isometric tension Piazzesi u. a. (2007); unloaded velocity of contraction must be comparable to experimental value. Moreover, in the MS model we also require that in the relaxed muscle about 5 % of myosin motors are in the ON state, as hypothesized in Linari u. a. (2015), and that the mechanosensing recruitment reaches the saturation at T greater than 0.5 T 0 . The force velocity curve simulated by the conventional model (Figure 3 A, red squares) shares the common limit described in the introduction: despite both the maximum velocity and the rising phase after activation (Figure 3  The limits of the conventional approach are clearly detectable when the force is allowed to recover isometrically after a period of unloaded contraction (Figure 3 B). The time constants, for the 20 ms (green line) and the 40 ms (red line) drop in tension, are equal to each other, as expected by a model without a feedback from the tension. With the current parameters, the simulated time constant is comparable to the one observed for the 20 ms drop but its pretty higher than that of the 40 ms drop. A different choice of parameter would modify the outcome, but the two rates can be made different only introducing further hypotheses in the conventional model.
The MS model includes the stress dependent transition between the two detached states. It is able to fit excellently the force-velocity curve (Figure 3 A, blue dots), while keeping the proper fitting for the rising phase during isometric activation (Figure 3 C, blue curve). Consequently, the power output at intermediate level of velocity of contraction is preserved at the physiological values. The model is able to keep a high maximum power thanks to a relatively high attachment and detachment rates, while the rising phase is slowed by the ON-OFF transition. These results strongly support the idea that the mechanosensing mechanism in the recruitment of ON or active myosins may solve the long standing conflict described before.
The simulated tension transients after unloaded contraction are shown in Figure 3 C. The curves for the tension development after initial activation, and after an imposed shortening at zero tension for 20 and 40 milliseconds, closely follow the single exponential fitting with time constants respectively of 34 ms, 24 ms and 28 ms. These values are the same as the ones observed experimentally Linari u. a. (2015), and also  2007)), conventional model: red squares, MS model: blue dots. B Conventional model. Simulated rise in tension after initial activation (blue), and tension recovery after 40 (red) and 20 (green) milliseconds at zero tension shortening phase. The rise in tension is comparable with the experimental data, as well as the recovery after 20 ms at zero tension, but the recovery rate does not depend on the duration of the zero tension period. C MS model. Simulated rise in tension after initial activation (blue), and tension recovery after 40 (red) and 20 (green) milliseconds at zero tension shortening phase. The single exponential fitting with the experimentally observed time scales well match the simulated tensions (τ R = 34 ms for the rise after activation, τ 10 = 28 ms after 10% of muscle length shortening, and τ 5 = 24 ms after 5% of muscle length shortening). Also the time delay are quite similar to the experimental data Linari u. a. (2015) (18 ms, 8 ms, 1 ms for the activation, 40 ms and 20 ms protocols respectively. Experimental data are 21 ms, 2 ms and 0 ms). the time delays are quite similar (see caption of Figure 3. Interestingly, the characteristic timing of the observed mechanosensing mechanism, is compatible with the attachment and detachment rate constants required to fit the high power output at intermediate velocities. Differently to the conventional model, now we can reproduce a lower time constant after longer zero-tension periods, because an higher number of myosin motors move to the OFF state in MS model, thanks to the mechanosensing system. The ON-OFF transition can be satisfactorily simulated as model mimics the increase of the myosins in OFF state during ramp shortening, roughly following what observed experimentally ( figure 4). It is important also to test the predictive ability of the model in relation to another behavior which has already been connected to the above discussed conflict, the so-called fast recovery of the power stroke. This behavior has been shown first about 25 years ago by Lombardi and co-workers Lombardi u. a. (1992). The tension recovery after a small "test" step in length, which follows a bigger "conditioning" step during steady state of the isometric contraction, increases with the time interval between the two steps, showing a greater recovery in tension already after 4 ms. The recovery increases more and more up to 15 ms after the first step. This evidence is incompatible with the slower time scale of the attachment detachment process, as inferred during the rising phase in conventional models. In Figure 5 are shown the simulations of the experimental protocol used for the multiple steps in Lombardi u. a. (1992). The model is able to almost quantitatively fit the experimental behavior, clearly showing an improved capacity to recover the original tension with the delay of few milliseconds from the conditioning step.
Finally, the analysis of the energetic predicted by the model requires some special considerations. As described before, in both the conventional and the MS models, the motors can detach from any attached state, also in the pre-power stroke. This event may be prevented in real muscle, or it may be possible but without the associated ATP splitting. Moreover, while the F atp drop imposed in E c (16 KT ) is lower than total ATP energy ( about 20 KT ), the H57 hypothesis implies an energy consumption which is difficult to asses. Both effects affect the efficiency predicted by the model. Anyway, associating one ATP splitting for each detachment event, the simulated energetic is compatible with data in the literature (Figure 6).

Conclusions
In this study we have shown that the tension dependent recruitment of active motors from the OFF state makes the tension development after activation less sensitive to the acto-myosin cycle rate constants and thus allows higher rate constant of attachment. This mechanism, which is substantiated by recent experimental findings Linari u. a. (2015), solves the long standing discrepancy between the high power output and the slow rise in tension in skeletal muscle contraction modeling. This conflict, surviving since the first model based on molecular mechanism proposed by A.F. Huxley on 1957HUXLEY (1957, has been previously solved only introducing ad-hoc hypotheses on the actin or myosin properties. We have also shown that the experimentally observed kinetics of this mechanism are quantitatively in agreement with our hypothesis. Perturbing in some way ON-OFF kinetics would allow to experimentally test this hypothesis. For instance, increasing, chemically or mechanically, the number of active motors in the relaxed muscle, would generate a faster tension rise after activation without affecting the force velocity relationship. Alternatively, perturbing the tension in the thick filament during activation would modify the number of active motors, and consequently affect the mechanical response of the muscle. The ON-OFF regulation, requires a tension feedback which is likely to be generated by a molecular component other than myosin or actin, like titin or the C-protein. How this feedback information is given to each myosin motor has to be explored experimentally. The tension may modify the state of the non-myosin proteins in the thick filament, affecting the rates of each motor in the same way. Or the extension of the thick filament may disturb the myosin-myosin interaction in the OFF state, inducing a nearest-neighbor cooperative interaction. The latter is required to describe the steep slope in the pCa-force curve Rice u. a. (2003) and the mehcanosensing mechanism may then contribute to the known coupling through the tropomyosin filament. Further experimental evidence are required for an explanation of the tension dependence, which is imposed on phenomenological basis in this paper. Nevertheless, in this work we have shown that including a new, experimentally well supported, structural property of the thick filament in muscle modeling is as important as the precise characterization of actomyosin interaction.

Supplementary Information 1 Supplementary text
A Monte-Carlo approach is used in both models, where the state and position of each myosin is followed at each time step ∆t = 10 −6 s. At each time step, the algorithm updates the position x i of the i-th myosin with respect to its anchor on the thick filament, considering the current position of the actin filament, the coupling factor for the myosin motors. In this way, the total tension generated is given by: where N XB is the total number of myosin heads, and N f il is the number of filaments in the sarcomere, while k is the stiffness of the myosin motor. The position of the myosin head is first updated based on the position of the actin filament, z, in the previous time step: Here z a i and x a i are the positions of the actin filament and the myosin pre-stretch, respectively, at the time of attachment of the i − th myosin head. n ps equals 0, 1 and 2 when the myosin head resides in the S0, S1 and S2 minimum, respectively. Having updated the position of each myosin head, the algorithm updates the position of the actin filament by the Newton-Raphson method. In each iteration, the total force is computed by Equation S1.
In both models, rate constants in the attached state are computed through a Kramers-Smoluchovsky approximation of the corresponding Langevin equation of a diffusing, over-damped particle with drag coefficient η in the sinusoidal energy landscape E c (x) + 1/2k(x − x 0 ) 2 (see Matherials and Methods): ω is zero when the myosin head is detached and one when it is attached. Γ(t) is a random term satisfying < Γ(t) >= 0 and < Γ(t 1 ), Γ(t 2 ) >= 2δ(t 1 − t 2 ) (white noise).
The algorithm for the model requires a pre-step outside the time loop:: 0. Generate a matrix of rate constants for each stretching level between -80 nm and 20 nm through a numerical integration of the Kramers-Smoluchivsky approximation of the Langevin equation describing the motion of a particle in the defined potential energy. Then: 1. Compute the tension acting on each thin filament and computing the new value of k 01 2. Define the new stable state of each myosin (OFF, ON, Pre or post power stroke) using the rate constants for the transitions k 01 , k 10 , k 12 , k 21 ,k f and k b .
3. Update the positions of each myosin head and the actin filament, depending on the setup and external conditions, and calculate the total force. 4. Update the position of the actin filament through an implicit method 5. Increment the time by one time step and return to Step 1 until total time is less than the prescribed