Theoretical efficiency limits and speed-efficiency trade-off in myosin motors

Muscle myosin is a non-processive molecular motor that generates mechanical work when cooperating in large ensembles. During its cyle, each individual motor keeps attaching and detaching from the actin filament. The random nature of attachment and detachment inevitably leads to losses and imposes theoretical limits on the energetic efficiency. Here, we numerically determine the theoretical efficiency limit of a classical myosin model with a given number of mechano-chemical states. All parameters that are not bounded by physical limits (like rate limiting steps) are determined by numerical efficiency optimization. We show that the efficiency is limited by the number of states, the stiffness and the rate-limiting kinetic steps. There is a trade-off between speed and efficiency. Slow motors are optimal when most of the available free energy is allocated to the working stroke and the stiffness of their elastic element is high. Fast motors, on the other hand, work better with a lower and asymmetric stiffness and allocate a larger fraction of free energy to the release of ADP. Overall, many features found in myosins coincide with the findings from the model optimization: there are at least 3 bound states, the largest part of the working stroke takes place during the first transition, the ADP affinity is adapted differently in slow and fast myosins and there is an asymmetry in elastic elements.

Introduction we only constrain a small number of parameters that are physically limited (like the stiffness and the nucleotide binding rates) and determine all other parameters by numerical optimization.

Model and methods
In our description of the non-processive motor we follow T.L. Hill's formalism [37], which is the basis of most myosin models. The duty cycle consists of attaching to the track (actin filament), a conformational change in the lever (working stroke), release of hydrolysis products and binding of a new ATP molecule followed by detachment, reverse conformational change (recovery stroke) and ATP hydrolysis. We describe the cycle of the motor with N S chemical states, of which N B are bound states (Fig 1A). Each state represents a unique conformation, characterized by the unstrained lever position d i and the free energy G i . With X we denote the position of the track relative to the backbone (thick filament), with x A the position of the binding site on the track and with x M the unstrained position of a motor on the backbone (Fig 1B). The strain on a motor is X + x A − (x M + d i ) = x − d i and its elastic energy U i (x) = U(x − d i ). In the simplest case of a harmonic potential, U(x) = Kx 2 /2 with a spring constant K (Fig 1C). The transitions are stochastic with rates

PLOS COMPUTATIONAL BIOLOGY
Theoretical efficiency limits and speed-efficiency trade-off in myosin motors The reverse transitions have the rates For transitions that involve attachment (state i unbound and state i + 1 bound) we set α i = 1, meaning that the strain only affects the attachment rate, while the detachment rate is constant. The strain dependence of the attachment rate is indirect and reflects the fact that in order to reach a binding site, the elastic element has to stretch through thermal fluctuations. Because of the distances involved, we consider this dependence stronger than the strain dependence of the detachment rate and neglect the latter. Likewise, we set α i = 0 for detachment transitions. If a transition involves binding of ATP, ADP, or Pi, the corresponding rate k 0 is k þATP ½ATP�, k þADP ½ADP� or k þPi ½Pi�, respectively.
For each state i we define P i (x) as the probability density to find the motor in that state with strain x. The probabilities are normalized to When the filaments move with velocity v = −dX/dt against each other, the strain on the motors changes according to dx/dt = − v. The system is described by the following set of master equations with the probability flux densities In detached states, the positions are quickly thermally equilibrated, and the corresponding master equations are @ @t P T i ðtÞ ¼ State i + N S is identical to state i, but because one ATP molecule has been hydrolyzed during the cycle, its free energy is G iþN S ¼ G i þ DG ATP . For a constant velocity v, the probability densities are determined by the stationary solution of Eqs (3) and (6). Then the total flux is conserved around the cycle and corresponds to the total rate of ATP hydrolysis per motor, The average force per motor follows as where U 0 i ¼ @U i =@x is the spatial derivative. To understand the origin of dissipation in the working cycle of myosin, we can compute the entropy production rates for individual transitions. At position x, each state can be assigned a chemical potential The transitions are generally out of equilibrium and lead to entropy production with density The total dissipation rate is obtained by summation over all transitions and integration over all positions x. One can verify that the total dissipation (time derivative of the total entropy in the system and the heat bath) equals the free energy of ATP hydrolysis per unit time, reduced by the power output [38]: In the second line we re-ordered the terms in the sum, in the third we applied Eqs (3) and (9) and in the last we performed partial integration. The efficiency follows as The model equations can be non-dimensionalized by expressing all energies with the thermal energy k B T and all distances with the power stroke size d. The velocity is expressed with k max þATP ½ATP�d, where k max þATP ½ATP� is the rate of ATP binding. While most model parameters are obtained by numerical optimization, there are three dimensionless parameters that are physically restricted. First is the dimensionless free energy of ATP hydrolysis À DG ATP =k B T, which is close to 25 under physiological conditions [5]. The second is the dimensionless stiffness which is constrained because of the elastic nature of a protein structure. We express it as the elastic energy of the spring when displaced by the working stroke distance d, in units of thermal energy: 1 2 Kd 2 =k B T. For myosin we estimate K ¼ 3 pN=nm and d ¼ 8 nm [39], which gives 1 2 Kd 2 =k B T ¼ 24. In addition, in models with three bound states (N B = 3) and a prescribed velocity, the binding rate of ADP is also a fixed parameter and expressed with the dimensionless rate k þADP ½ADP�=ðk þATP ½ATP�Þ. Furthermore, the efficiency is independent of the motors' duty ratio (fraction of time it spends in the bound state), because the dwell time in the detached state slows down the cycle and proportionally reduces the force and the ATP consumption, but does not affect their ratio. Therefore, the relative free energy level of the detached state, G 0 − G 1 does not appear as a relevant optimization parameter. In cases where we only discuss the maximum efficiency, independent of velocity, the dimensionless velocity also becomes an optimization parameter.

Numerical solution
For a given set of parameters we determined the stationary solution (@P i /@t = 0) of the system of coupled linear differential equations given by Eq (3) with the boundary condition P i (x ! 1) = 0 for i = 1. . . N B . We carried out the numerical integration in negative direction in x using a non-adaptive 3-step backward differentiation formula (BDF) method, which is suited for handling stiff ODEs. Simultaneously, we integrated the rate of ATP consumption given by Eq (7) and force by Eq (8) (see S1 Source Files). The resulting efficiency is calculated as Z ¼ vF=ðÀ r ATPase DG ATP Þ. The efficiency as a function of the parameters that are listed in Table 1 for each case is used by a numerical optimizer using a quasi-Newton algorithm (routine E04JYF from the NAG Library, Numerical Algorithms Group). Solutions that only considered operating points outside a hysteresis in the force-velocity relation were obtained using a sequential quadratic programming procedure with nonlinear constraints (NAG routine E04UCF). Parameter scans were made by initializing the optimizer with the previous solution, once in ascending and once in descending direction, and additional verifications with different initial parameter values were carried out. The values of optimization parameters for all solutions shown in this paper are given in S1 Data.

Validity of the stationary solution
Our numerical solution assumes that the motor ensemble moves at a constant velocity during contraction. In myosin, this stationarity has often been questioned. In fact, some studies show that under high loads (close to stall) myosin motors can synchronize their cycles, leading to a step-wise contraction [40][41][42]. Signatures of such steps can be seen in muscle transiently after 0 (no effect on efficiency)

PLOS COMPUTATIONAL BIOLOGY
Theoretical efficiency limits and speed-efficiency trade-off in myosin motors changes of load [43]. Kaya et al. reported step-wise motion in small groups of myosin motors [44] under high load. In these studies, the coordinated stepping appears in a narrow regime close to isometric conditions. We therefore do not expect it to affect the maximum efficiency, which is achieved at a lower load and higher velocity. Furthermore, recent direct observations with high-speed AFM (atomic force microscopy) showed no cooperativity between myosins on non-regulated thin filaments [45].
We verified the validity of the stationary solution by comparing our numerical solutions with stochastic simulations on a finite ensemble of motors. For a selection of optimal parameter sets, we simulated a group of N m motors pulling against a constant load. The simulations were carried out with a Gillespie algorithm and rapid mechanical equilibration at each step [42]. All results show a very good convergence for two-digit motor numbers (S1-S4 Figs). The number of motors needed is similar to the numbers that were needed to achieve a high efficiency in single-filament experiments [46]. No coordinated steps are visible in the simulated traces. The efficiencies obtained from the simulation converge towards the results of the stationary solution, and the velocity becomes constant in time for large motor numbers.

Velocity-independent maximum efficiency
We first numerically determine the maximum efficiency η max for a system with N B = 2, 3, . . . bound states. Because the velocity is unconstrained we always evaluate the efficiency at the peak of the η-v relationship (see Fig 1D for an example). If we put no constraint on the kinetic rates, η also depends on the elastic constant 1 2 Kd 2 =k B T, the free energy differences between bound states G i /k B T, the ratios between kinetic constants k i /k j , coefficients α i and for N B > 2 on the lever position in the intermediate states d i for 2 � i � N B − 1. We also fixed the parameters that do not influence the optimization outcome (i.e., one kinetic constant, a constant offset in free energy levels and the free energy difference between unbound and bound steps). All other parameters have been obtained by a multidimensional optimization procedure. All parameters of the model are listed in Table 1.
The results are shown in Fig 2. The dependence of maximum efficiency on the free energy of ATP hydrolysis is non-monotonic (Fig 2A). If the ATP, ADP and Pi concentrations are close to chemical equilibrium, −ΔG ATP /k B T � 1, the efficiency is η max = 0.162 for N B = 2. A higher efficiency is possible if the nucleotide concentrations are further away from equilibrium. A similar effect, namely an efficiency that is maximal outside the linear response regime, has already been observed in some scenarios with Brownian ratchets [47]. For −ΔG ATP /k B T = 25 (physiological value), the efficiency limit is 42%. However, this efficiency is only achievable if the motor ensemble is working with a prescribed velocity. Namely, the maximum efficiency is found in a region with an anomalous force-velocity relationship ( Fig  1D). In a more common fixed-force scenario, the ensemble cannot be held at that point without jumping also to a less efficient point with the same force. In the fixed-force regime, we therefore impose a constraint that the operational point has to be outside the hysteresis. This reduces the maximum possible efficiency to 38% (thin solid line in Fig 2A). The origins of dissipation are analyzed in Fig 2C, which shows the rate of entropy production (Eq (10)) for each transition and each head position x (equivalent of strain). The losses mainly result from a premature detachment of strained heads (red line in Fig 2C). The model with N B = 3 bound states, on the other hand, can already achieve very high efficiencies, but only with very stiff springs (dotted lines in Fig 2A). Here the attachment/detachment of heads occurs primarily from/to the position close to the unstrained state (Fig 2D). The maximum efficiency as a function of the elastic constant K is shown in Fig 2B. One expects that the elastic energy stored into the spring during the working stroke imposes an upper limit on the useful work, which gives an efficiency of Z � 1 2 Kd 2 =jDG ATP j. In the following, we provide an analytical argument for the efficiency bound in the limit of very soft springs, Kd 2 � |ΔG ATP | (S5 Fig). In this limit the forces and the elastic energies become small and lose the influence on the transition rates. The rates can be regarded as irreversible. Models with strain-independent transition rates can be solved analytically [31]. The unloaded velocity is maximal when the first transition is fast k 1 ! 1 and contains the full working stroke (d 1 , with the remaining transitions having equal rate con- 48]. The maximum force per motor is Kdt on /t cycle , while the ATPase rate is always 1/t cycle . The maximum power and efficiency are achieved at 1/2 the stall force and 1/2 the maximum velocity. At that point, the work per ATP is 1 2 Kd 2 ð1 À 1=N B Þ and the efficiency is 1/N B ). Therefore, the maximum efficiency in the limit of soft springs grows linearly with the stiffness with a proportionality constant that depends on N B (Fig 2B). The upper limit of 1 2 Kd 2 can only be reached with an infinite number of bound states. This solution reveals another interesting observation, namely that the transitions between bound states are not necessarily fast in the optimal regime. Although a slow transition increases dissipation, a cascade of transitions with similar rate constants narrows the distribution of times between attachment and detachment, i.e., it makes it more deterministic (S5(D) Fig). A narrower distribution means that the elastic energy of strained springs just before detachment, which is another source of dissipation, will be smaller. In the extreme limit of fully deterministic attachment times, it is even possible that all motors detach with exactly zero strain. This result is in contrast with processive motors, where increasing the kinetic constants always improves the efficiency at a given velocity [22].

Maximum efficiency at a given velocity
In the second scenario we determine the maximum efficiency at a given velocity v while imposing upper limits on two rate constants: the rate of ATP binding with detachment k max þATP and the rate of ADP binding k max þADP . Both rate constants are of a similar order of magnitude, * 1 μM −1 s −1 , in most myosins [49,50], but because of the lower ADP concentration [5], we expect a ratio of k þADP ½ADP�=k þATP ½ATP� � 1=100. The maximum efficiency for different ratios between the two limiting rates is shown in Fig 3A. It displays a trade-off between speed and efficiency. In the limit of a very slow sliding velocity the maximum efficiency with constrained kinetic rates reaches the values from the previous section, independent of the ratio between the bounds on kinetic rates. The fast transition between the first bound states contains most of the working stroke (d 2 � d 1 ). The optimal stiffness is at its upper limit which allows for the high exerted force and a better control on the strain dependent detachment rate. For high sliding velocities the maximum efficiency is decreased and a large part of the working stroke is allocated to the transition with the limited kinetic rate (d 2 * d 1 ) The optimal stiffness becomes softer to reduce negative forces caused by the imprecise position of attachment and detachment of heads.
An interesting outcome of the optimization is that the free energy difference between the last two bound states G N B À 1 À G N B depends strongly on the velocity for which the motor is optimized. At low velocities, an uphill transition slows down the detachment and augments its strain sensitivity, see the potentials G i + U i (x) in the top panel of Fig 2D. At the same time, a low energy of the penultimate bound state increases the energy difference available for the working stroke. On the other hand, motors optimized for high velocities have a strong downhill transition before detachment, see the top panel of Fig 3C. In this case, quick detachment of motors is more important than the precise control over the position where the detachment takes place. In myosin the last transition between the bound states is the release of ADP. The free energy difference of this transition can be expressed with the ADP dissociation constant The resulting optimal K ADP D as a function of the velocity is shown in Fig 3B. It is noteworthy that the dependence is similar in both cases, for rate limiting ATP-(black curve) and ADP-binding (other curves). The finding is in agreement with the experimental observation that ADP affinity is the main mechanism of adaptation between fast and slow muscles [49]. Over a velocity range of 0.3 to 7 μm/s, the K ADP D scales approximately with *v 2 [51]. The predicted optimal ADP affinities are remarkably close to the quadratic law (dashed line in Fig 3B) over a wide interval of velocities.

Anharmonic elastic potential
In the previous section, we saw that the stiffness of the elastic element had two conflicting effects. For motors at the beginning of the power stroke, a stiffer motor can produce more force. In a stiffer motor, transitions also take place in narrower intervals, further reducing the dissipation. On the other hand, motors that remain bound past x = 0 induce an effective drag. This raises the question whether efficiency can be improved by allowing an anharmonic spring. Asymmetric models for the myosin elasticity have been proposed previously in order to explain the measured properties of myosins [44,[52][53][54] or directly observed [55]. We have therefore relaxed the condition that the potential U(x) be harmonic and allowed for an arbitrary shape, provided that 0 � U"(x) � K. Numerically, we have divided the potential into a finite number of segments and run the optimization with the stiffness in each segment as a separate (constrained) parameter. The results show, however, that in all cases analyzed the optimal stiffness has a bimodal shape, reaching the constraint at x > x a and being zero (U" = 0) for x < x a , where x a is a transition point obtained by optimization.
The dashed lines in Fig 3A show the maximum efficiency if we allow anharmonic potentials. We observe an interesting transition: for slow velocities, a harmonic potential with the highest allowed stiffness is still optimal. However, at higher velocities, the optimal potential is asymmetric (Fig 3D). The shape of the potential reduces the negative forces caused by the post-powerstroke heads once they pass the potential minimum. At the same time, the remaining barrier, which is of the order of thermal energy, keeps the detached heads close to their unstrained position. In muscle myosin, such a potential can result from the elastic buckling of the tail domain and has indeed been observed in single molecule experiments [39,56]. In finite ensembles of motors, asymmetric stiffness can have the consequence that the velocity becomes limited by the attachment, rather than the detachment rate [54,57]. The softening of the spring in motors optimized for high velocities, on the other hand, is not seen in myosins. In fact, several comparisons between muscle types showed a higher stiffness in faster isoforms [58][59][60][61]. We note, however, that the measured stiffness is largely determined by the pulling motors (x > 0) and drawing conclusions on the regime with negative strain (x < 0) is difficult.

Conclusions
We looked into the problem of myosin energetics from the reverse perspective: how would an optimally designed myosin motor work? We only consider a few physical constraints that cannot be arbitrarily altered in the course of optimization. These include the number of bound states, which is largely tied to the ATP hydrolysis cycle. We further constrain the maximum stiffness, which we expect to be limited by the elastic properties of a protein. Finally, the second order rate constants of ATP and ADP binding are broadly conserved between myosins, which suggests that they are also close to their physical limits. Besides those constraints, our model makes some further assumptions. For one, we assume a single working cycle without any side branches, for example detachment in the ADP state. Such side branches generally take place between states with a higher free energy difference and therefore increase the dissipation (they can also be seen as leaks in the cycle). On the other hand, we neglect the discrete nature of actin subunits and assume that the myosin heads can bind anywhere on the actin filament. While the effect of discrete binding sites has been estimated as small [42], we expect that restricting the sites will somewhat increase the losses of the transition involving initial binding (0 ! 1). The assumptions that the detachment rates are independent of force and that the other transitions depend in the simple way on the free energy difference present a restriction in the space of possible models. Furthermore, like virtually all myosin models, we assume that the stiffness of the elastic element is the same in all states. Also, we do not consider other optimization criteria that may be physiologically relevant, such as the maximum (isometric) force per motor or the ability of motors to work efficiently over a broader range of velocities. Beyond these constraints and assumptions, all parameters are the result of optimization. We identify three limiting factors for the efficiency: the number of bound states, the stiffness and the kinetics of ATP and ADP binding. Efficient motors require at least 3 bound states. The optimal free energy difference between the last two bound states (in myosin: ADP and without a nucleotide) depends strongly on the velocity for which the motor is optimized. Slow motors require an uphill transition, while fast ones work better with downhill. Another adaptation concerns the stiffness: at low velocities, a stiff potential allows them to exert a high force and to better control the strain dependent detachment rate. The importance of stiffness for a high efficiency was already highlighted in several muscle models [40,41,62,63]. At high velocities the optimal potential becomes asymmetric: stiff for pulling forces and compliant in the pushing direction, above a certain threshold force. Many of these features are found in myosins: the chemical cycle has 3 distinct bound states [64], the largest power stroke takes place at the beginning of the cycle (connected to the release of Pi [65]), the energetics of the ADP release is the main source of variability between slow and fast myosins [51] and the elastic properties of the tail lead to an asymmetric potential [39]. A cross-comparison between different muscle types also shows a trade-off between speed and efficiency [51], as expected from our calculations. Similar observations were also made about a trade-off between power density and efficiency [66]because the power output depends on more parameters than just the speed, it is more difficult to compare. Using concepts from stochastic thermodynamics, we were also able to identify the sources of dissipation inside the cycle, which can be more informative than just the efficiency of the cycle as a whole. We show that the largest free energy losses take place in the steps related to attachment and detachment from the actin filament. This shows that the rules governing the efficiency of non-processive motor proteins are very different from the processive ones like kinesin or F1-F0 ATP synthase, which were frequently studied from the perspective of stochastic thermodynamics [13,22,24,25,67,68], and that the stochastic nature of attachment and detachment events leads to unavoidable losses.

PLOS COMPUTATIONAL BIOLOGY
Theoretical efficiency limits and speed-efficiency trade-off in myosin motors