Fluctuating Currents in Stochastic Thermodynamics II. Energy Conversion and Nonequilibrium Response in Kinesin Models

Unlike macroscopic engines, the molecular machinery of living cells is strongly affected by fluctuations. Stochastic Thermodynamics uses Markovian jump processes to model the random transitions between the chemical and configurational states of these biological macromolecules. A recently developed theoretical framework [Wachtel, Vollmer, Altaner:"Fluctuating Currents in Stochastic Thermodynamics I. Gauge Invariance of Asymptotic Statistics"] provides a simple algorithm for the determination of macroscopic currents and correlation integrals of arbitrary fluctuating currents. Here, we use it to discuss energy conversion and nonequilibrium response in different models for the molecular motor kinesin. Methodologically, our results demonstrate the effectiveness of the algorithm in dealing with parameter-dependent stochastic models. For the concrete biophysical problem our results reveal two interesting features in experimentally accessible parameter regions: The validity of a non-equilibrium Green--Kubo relation at mechanical stalling as well as negative differential mobility for superstalling forces.


I. INTRODUCTION
Understanding the complex biochemical processes which are responsible for cellular metabolism is one of the key questions in modern biophysics. The quantitative analysis of so-called molecular motors, which are the small machines transforming different forms of energy into one another, is at the center of these efforts [1][2][3]. In recent years scientists developed techniques that allow the systematic observation and manipulation of these biological macromolecules [4]. Under in vivo conditions, (electro-)chemical gradients in the cell maintain these systems out of equilibrium. From a thermodynamic perspective, one is interested in the currents of heat, matter and energy that flow through a molecular motor, because they allow, for instance, the definition of its efficiency.
In analogy to macroscopic engines, molecular motors are described by thermodynamic cycles in a space of biochemical and configurational states. In contrast, the energy scales involved in biochemical energy conversion are only a couple of times larger than the thermal energy. Consequently, thermal fluctuations cannot be neglected, and their dynamics must be modeled as a stochastic process that reproduces the stochastic time series observed in experiments. Stochastic Thermodynamics refers to a general framework for a consistent definition of fluctuating work and heat currents on the level of the fluctuating time-series [5,6]. The common model for molecular motors are dynamically reversible Markov jump processes, which can be thought of as (memoryless) random walks on a biochemical network of states [3,[7][8][9]. In an accompanying publication [10] we investigated the asymptotic statistics of such systems from the perspective of the cycle topology of the network of states. In particular, we developed an efficient method to calculate all cumulants of arbitrary fluctuating currents analytically.
Here, we are interested in the first and second order fluctuation statistics, i. e. the expressions for macroscopic average currents (like the motor's velocity) and Green-Kubo time-correlation integrals (like its diffusion constant). To be concrete, we use the analytic nature of our method to analyze the parameter space of different stochastic models for the motor protein kinesin [11][12][13], which were designed to reflect typical force-spectroscopy experiments [14][15][16][17]. Besides illustrating the insights that thermodynamic cycles provide into the motor dynamics, our results uncover interesting model predictions and thus indicate directions for future experimental research: The validity of a nonequilibrium fluctuation dissipation relation at mechanical stalling as well as negative differential mobility, commonly referred to as "getting more from pushing less" [18].
This work is structured as follows. In Sec. II we briefly review the results of Ref. [10]. In contrast to the formal exposition there, here we focus on the implementation of a universally applicable algorithm for the efficient calculation of averages and correlation integrals of fluctuating currents in Stochastic Thermodynamics. Sec. III thoroughly discusses how to apply our universal method in the concrete biophysical context of a kinesin model. In Sec. IV we give detailed account of kinesin's chemical (ATP hydrolysis) and mechanical (displacement) currents as functions of their conjugate chemical and mechanical drivings. We conclude in Sec. V with a discussion of the main conceptional and biophysical insights.

II. FLUCTUATING CURRENTS AND THEIR STATISTICS
In this section we introduce our mathematical notation and -based on the general results presented in Ref. [10] -provide a concise recipe for calculating the averages and asymptotic (co-)variances of two fluctuating currents in a dynamically-reversible Markov process on a finite state space. Such averages and covariances play a major role in Stochastic Thermodynamics [8,19,20]: They correspond to physical steady-state currents and timecorrelation (Green-Kubo) integrals [21]. To be concrete, we exemplify topological concepts for both a four-state and a six-state Markov process, Fig. 1. In Sec. III we interpret these examples as models for the molecular motor kinesin.

A. Currents for Markovian processes
Stochastic processes on a finite state space V = {v 1 , v 2 , . . . , v N } are called Markovian jump processes. Henceforth, we consider the time-continuous and homogeneous case. A realization, or trajectory, of the process starting at time t 0 = 0 in a state γ 0 ∈ V is a collection of jump times t k > 0 and visited states γ k ∈ V with k ∈ N. Transitions from a state v i ∈ V to a different state v j ∈ V occur at a given constant rate w i j . For thermodynamic consistency [6-8, 20, 21], we require dynamical reversibility, i. e. w i j > 0 ⇔ w j i > 0 . With this constraint we draw the state space as an undirected graph G with the N states as vertices and M admissible transitions as edges, cf. Fig. 1. In addition to dynamic reversibility we assume that the state space is connected, which ensures ergodicity of the process.
An ensemble of trajectories with initial probability distribution p(0) = (p 1 (0), . . . , p N (0)) on V evolves according to the master equation [22]: d dt p(t) = p(t)W or, in components, where we use the convention w i i = − j =i w i j . Ergodicity of the process implies that there is a unique steadystate probability distribution π satisfying 0 = πW to which all initial conditions will converge eventually. The quantities J i j = π i w i j − π j w j i represent the steady-state probability currents between two states v i and v j . The special condition where all the steady-state currents vanish, J i j = 0, is called detailed balance or equilibrium. Since we are interested in the currents of non-equilibrium systems, we will not assume detailed balance in the following.
In order to account for general currents, e. g. changes in energy, entropy, particle numbers or physical position, we introduce jump observables. A jump observable ϕ assigns a weight ϕ i j to the transition from v i to v j , where we require anti-symmetry: ϕ i j = −ϕ j i . The macroscopic current associated to a jump observable ϕ is In order to illustrate the concept we provide two examples. For a pair of states (v i , v j ) we define a simple but important case of jump observable: the counting observable ϕ (i,j) counts +1 for a transition from state v i to v j and −1 for a reverse transition from v j to v i . To every other transition it associates a weight of 0. In this case the macroscopic counting rate from v i to v j equals the steady-state probability current between these states: c(ϕ (i,j) ) = J i j . This expression obviously vanishes if transitions between v i and v j are impossible. In Ref. [10] we emphasized that counting observables form a basis of the space of jump observables: every jump observable can be expressed as an appropriate linear combination of counting observables.
Another important example is the dissipation in stochastic thermodynamics. It is derived from the jump observable σ that takes the values σ i j = ln is non-negative and vanishes only at equilibrium, i. e. if and only if we have detailed balance. In Ref. [8], Schnakenberg discussed an analogy between Markov jump processes and Kirchhoff's laws for electric circuits. Kirchhoff's current law states that the steadystate currents balance at each vertex v i , i J i j = 0 . The fundamental object of Schnakenberg's network theory is a set of B = M −N +1 fundamental cycles {ζ }. They are obtained from a spanning tree for the graph (cf. Fig. 1 as well as Sec. II B). For now, think of a fundamental cycle ζ as a tuple of consecutive states γ 0 , γ 1 , . . . , γ m( ) = γ 0 that form a self-avoiding and closed trajectory. Cycles are defined up to cyclic permutations. Adding the contributions of a jump observable ϕ along the transitions in a fundamental cycle ζ gives its circulation Then Kirchhoff's voltage law states that the circulation for the voltage drops in any electrical circuit vanishes. Another result obtained in this context is the Schnakenberg decomposition of the macroscopic dissipation rate [8]: where c is the steady-state probability current associated to the cycle ζ [10]. The circulationsσ of the dissipation are called the cycle affinities. In the context of irreversible thermodynamics [23], the Schnakenberg decomposition (3) identifies the cycle affinities as the generalized forces which are conjugate to the cycle currents c .
In an earlier publication [13] the authors pointed out that Schnakenberg's decomposition is equally applicable to other observables ϕ. As a corollary, one may express the average current c(ϕ) = c φ using only the cyclic structure of the graph. Fig. 1 shows two graphs with four and six states, respectively, but having the same cyclic structure. Collecting all B circulations of an observable ϕ in a vector gives its chord representation ϕ H := (φ 1 , . . . ,φ B ) ∈ R B . The detailed mathematical background of this representation is discussed in Ref. [10].
Here we are interested not only in the macroscopic expectations of currents, but also in their higher order statistics. Consequently, the object of study in this work are fluctuating currents. The instantaneous current  ϕ (t) derived from a jump observable ϕ along a trajectory (t k , γ k ) is defined as thus accounts for the total change of the observable ϕ along a random trajectory with a random number n(T ) of jumps up to time T . Hence, this time-integral is a random variable with its own statistics. A typical realization of the time-integrated current, and thus its expectation value, grow linearly in time. Due to ergodicity, the time-averaged current ϕ T := 1 T ϕ(T ) converges to the macroscopic current c(ϕ) in the long time limit: Equivalently, one can average the fluctuating current over trajectories in the steady-state ensemble:  ϕ (t) =  ϕ (0) = c(ϕ) , where one can exploit the fact that the steady state is time independent. This ensures that the macroscopic current c(ϕ) is the expectation value, or mean, of the fluctuating current  ϕ (t).
Another important statistical measure are the correlations of two currents  ϕ and  ψ . A measure for this correlation is the Green-Kubo integral : Similar to the case of the average currents, ergodicity allows us to replace the steady-state ensemble average by a time average over the argument of the first current  ϕ . As a consequence [21], the correlation integral corresponds to a properly scaled covariance of the time-averaged currents: As such, the macroscopic current and the Green-Kubo integral are the first two scaled cumulants of the pair ϕ T , ψ T of time-averaged currents. Scaled cumulants are defined as derivatives of the scaled cumulantgenerating function, which can be obtained using methods from Large Deviation theory [10,21,24,25]. Higher order derivatives represent higher orders of the statistics, such as skewness and kurtosis.
In our accompanying publication Ref. [10] we prove a gauge invariance of the fluctuation statistics and show in detail how Schnakenberg's decomposition is extended to all cumulants of arbitrary observables. In the next section we give a brief review in form of a concise and efficient recipe for the first two orders.

B. Determining averages and (co-)variances of currents
The only ingredients needed for the calculation of the scaled cumulants are the transition matrix W and the jump observables representing the currents of interest. The elements w i j of the transition matrix as well as the jump observables may depend on the (physical) control parameters of the Markov processes in an arbitrary way. In the following, we present a simple yet efficient algorithm for the calculation of the first two scaled cumulants of two jump observables ϕ and ψ. It consists of three subalgorithms: topological, algebraic and physical. The first two steps are universal. Only the last step involves the jump observables in question. Note that at no step we explicitly calculate the steady-state distribution π, or the scaled-cumulant generating function. In most cases the latter is difficult or even impossible to obtain analytically.

Topology: defining fundamental cycles
The first step in the analysis addresses the topology of the graph G representing a network of states, cf. c. Identify the fundamental cycles: for every chord η ∈ H, its terminus and origin are connected by a unique directed path through the spanning tree. Adding the chord itself as a closure of this path results in the fundamental cycle ζ (Figs. 1b,d).

Algebra: determining the fundamental current cumulants
The second part of the algorithm involves the determination of the first two (joint) scaled cumulants of the fluctuating currents associated to the chords η ∈ H.
a. Write down the characteristic polynomial b. Identify the coefficients a 0 (q), a 1 (q) and a 2 (q) of χ H (λ; q) = N k=0 a k (q) λ k , i. e. the coefficients of the constant, the linear and the quadratic term.
c. Calculate the vector c =∈ R B with entries c = c(η ) and the scaled co-variance matrix C ∈ R B×B with entries C m = c(η , η m ), as follows (cf. Eqs. (12) of Ref. [10] for the derivation): where the partial derivatives ∂ a k := ∂a k (q) ∂q q=0 and the coefficients a k are evaluated at q = 0.

Physics: cumulants of jump observables
The concluding step of the algorithm finally characterizes the statistics of the fluctuating currents associated to the jump observables ϕ and ψ.
a. Sum the jump observables ϕ and ψ along the edges of the fundamental cycle ζ to obtain the circulations ϕ andψ . They are the coordinates of the chord b. The steady-state average of ϕ, and of the scaled covariance of ϕ and ψ then read We conclude the section with final remarks on the choice of the spanning tree in step 1, which is a priori arbitrary. Different choices yield different chords -and thus the different expressions for the fundamental current vector c and the fundamental covariance matrix C. However, expressions (9) and (10) are universal. In order to calculate the cumulants of any jump observable, no equations need to be solved. Via Eq. (9), any combinatorial complexity is hidden in (the derivatives of) the coefficients a k of the characteristic polynomial. The latter are calculated in a straightforward way either manually or by using a computer-algebra system. However, the final expressions (10) have fewer terms if some of the circulationsφ orψ along fundamental cycles vanish. It is thus worthwhile to take a careful look at the particular set of jump observables ϕ and ψ under consideration and choose a spanning tree that is optimal in that regard.

A. Kinesin as a molecular motor
Kinesin is a molecular motor which facilitates transport in eukaryotic cells. It moves along intracellular filaments called microtubules and plays a major role in many biological processes, including mitosis, meiosis and transport of cellular cargo. The most well studied variety of kinesin -both experimentally (see e. g. [4,[15][16][17] and References therein) and theoretically [11,12,[26][27][28][29] -is a protein dimer consisting of two identical subunits. Figure 2a shows a sketch of kinesin binding its intracellular cargo at its tail end. Kinesin's head end consists of two active sites which bind and unbind to the microtubule in alternating succession, thereby allowing the motor to perform mechanical steps of length L = 8 nm [30,31]. Due to the polarity of the microtubule, this motion has a preferred "forward" direction.
The energy necessary for this active directed transport is provided by the hydrolysis of adenosine triphosphate (ATP) into adenosine diphosphate (ADP) and inorganic phosphate (P). Unlike macroscopic motors, small molecular machines operate at low Reynolds numbers and inertia plays no role: chemical energy is not converted into mechanical energy by a transfer of momentum. Instead, kinesin's mechanical displacement is the result of a complex interplay of the strength of the microtubule binding at the active sites, which depends on their chemical composition. ATP-laden and empty sites bind strongly, while ADP-laden sites bind weaker [16,17]. Under physiological conditions the mechano-chemical interaction can be described by the "forward cycle" depicted in Fig. 2b [11]. Models that only treat the forward cycle feature tight coupling between the hydrolysis reaction and the stepping: each hydrolysis of an ATP molecule gives rise to exactly one motor step [14].

B. Experiments and models
An important biophysical question regards the force that kinesin generates for different concentrations of the chemicals ATP, ADP and P involved in the hydrolysis reaction. Typically, experiments measure this force by linking kinesin to a di-electric colloidal bead which resides in an optical trap [14][15][16]30]. Involved experimental setups allow the precise control of the pulling force F that the optical trap exerts on the motor against its typical direction of motion. The independent driving parameters are the non-dimensionalized force f := (LF )/(k B T ) and the non-dimensionalized chemical potential differ- Many experiments probe the stalling force f stall (∆µ), which is defined as the value of the force needed to bring the motor to a halt for a given chemical potential difference ∆µ. Under physiological chemical conditions, kinesin hydrolyzes ATP even at stalling forces [15,16]. The exact details of the kinesin stepping mechanism under high mechanical loads remain unknown and several models exist, cf. Refs. [11,12,17] and the references discussed in these publications. While these models differ in their details, they all feature more than only the tightly coupled forward cycle.
A prominent example of a thermodynamically consistent model was introduced in Ref. [11]. There, the key idea is to extend the forward cycle shown in Fig. 2b by the chemical compositions obtained from exchanging the trailing with the leading active sites, cf. Fig. 3a. In addition to the forward cycle F, the extended network features six states and has two additional cycles, Figs. 2b-d: the backward cycle B represents backward motion under hydrolysis, while the dissipative slip cycle D represents the futile hydrolysis of two ATP molecules without any stepping [27]. In such a multiple-cycle model, hydrolysis and mechanical displacement are not tightly coupled anymore. In Sec. IV A we will address the question of quasi-tight coupling, i. e. situations where the ratio of the average number of chemical and mechanical events predicted by the model is close to one.  Fig. 1c. Then, the fundamental cycles ζ1 ≡ F and ζ2 ≡ D are the forward and dissipative cycle (b) and (d), respectively. The backward cycle (c) is the linear combination B = ζ2 − ζ1, cf. Ref. [10].

C. Network theory for the kinesin model
In order to study quasi-tight coupling, energy conversion and the predicted response to changes in the driving parameters, we apply the algorithm presented in Sec. II B to the six-state model for kinesin, Fig. 3. For the first step of the algorithm, we choose the spanning tree and its chords in the same way as in Fig. 1c,d. Consequently, the fundamental cycles ζ 1 ≡ F and ζ 2 ≡ D correspond to the forward and dissipative cycles, respectively.
The second step of the algorithm requires the determination of the fundamental current vector c and the fundamental co-variance matrix C. With the enumeration of the vertices as in Fig. 3a the matrix W H (q 1 , q 2 ) reads: It is straightforward to write down its characteristic polynomial χ H (λ; q 1 , q 2 ) =: 6 k=0 a k (q 1 , q 2 )λ k , and to extract the coefficients a 0 (q) ≡ det W H (q 1 , q 2 ), a 1 and a 2 . Differentiating with respect to q 1 and q 2 , and evaluating at q 1 = q 2 = 0 yields the expressions ∂ a k appearing in Eqs. (9).
The third step requires the circulations of the jump observables of interest. For the present discussion, we consider the displacement d = ϕ (2,5) and the hydrolysis count h = ϕ (6,1) + ϕ (3,4) , which indicate a transition along the brown and purple edges in Fig. 3, respectively, cf. Sec. II A. Their matrix representations read The circulations of d and h simply count the number of the brown and purple edges in the fundamental cycles ζ 1 = F and ζ 2 = D, cf. Fig. 3b,d. Their chord representations thus are Note that the choice of the chords is optimal for the calculation of the present variables because one of the entries of d H vanishes (d 2 = 0), while this cannot be achieved for h H . After all, all cycles contain at least one hydrolysis event.
The corresponding macroscopic currents, i. e. the velocity c(d) and the hydrolysis rate c(h) are obtained from Eq. (10) as: Their scaled (co-)variances amount to In addition to displacement d and hydrolysis count h, we are interested in the jump observable σ i j = ln(w i j /w j i ) corresponding to the dissipation. As discussed in Sec. II A its circulations are the cycle affinities. The Hill-Schnakenberg conditions are necessary for the consistency of a Markov jump process with the thermodynamic notion of local equilibrium [3,9]. They state that the affinity of a cycle must express the (nondimensionalized) differences in the potentials of the reservoirs, cf. Refs. [7,8,19]. Upon completing the forward cycle F ≡ ζ 1 , an amount ∆µ of chemical energy is used by the system to perform a (dimensionless) amount −f of work against the pulling force. Similarly, a completion of the dissipative cycle D ≡ ζ 2 uses 2∆µ of chemical energy. Consequently, the chord representation of the dissipation reads σ H = (σ 1 ,σ 2 ) = (−f + ∆µ, 2∆µ).
The Schnakenberg decomposition thus lets us express the average steady-state dissipation by physical parameters and currents through fundamental chords D. The cycle perspective In the previous section, we expressed observable quantities only by means of their circulations around fundamental cycles and the fundamental first and second chord cumulants. Nowhere in these expressions do the number of states or the choice of a spanning tree appear explicitly. Hence, the same expressions are reproduced by any model with the same cycle topology -as long as the physics along the cycles, i. e.the circulations of observables, are the same. As exemplified in Sec II A in Fig. 1a,b, one can formulate a model on four states with the same cycle topology as the six-state model described in Fig. 3. Allowing only for single edges between states, such a four-state model is the minimal model featuring two independent cycles. An interesting question is how this model (and other reduced models) compare to more complicated ones. In Ref. [13] we used the idea of preserving the cycle affinities and circulations of jump observables along cycles (together with locality constraints) to develop a coarse-graining algorithm for stochastic models. Its application to the six-state kinesin model produced various topologically equivalent models, which all preserved the fluctuation statistics of the observables of interest almost perfectly. A disadvantage of this coarse-graining algorithm is that the individual transitions in the network of states lose their original interpretation.
In contrast, Fig. 4 shows a four-state model with a clear interpretation of the transitions. This has the advantage that the parameterization of the transition rates is found by the same physical arguments as the ones used in Ref. [11] for the six-state model. Details of the construction of this model are given in Appendix A. We will see in the next section that all the predictions of the six-state model are also found in topologically equivalent four-state models. This observation underlines the virtue of viewing models in ST from the perspective of cyclesan idea that was pioneered by Hill [7,19].

IV. RESULTS
One of the main messages of this work is that the algorithm presented in Sec. II B allows one to probe parametric models used in Stochastic Thermodynamics in a systematic way. In order to demonstrate the efficiency of our approach, we report on various non-trivial predictions of the six-state kinesin model. At the end of this section we will compare these results with other models.
Throughout this section, we choose the parameter range similar to Ref. [27] and vary −30 ≤ f, ∆µ ≤ 30. Then, in physical units the pulling force F varies between about −15 and 15 pN. Following Ref. [11], the chemical potential difference is adjusted by changing the ATP concentration while fixing the other chemical concentrations at physiological values (see also Appendix A). The physiologically relevant region for the chemical driving parameter is limited to about 20 < ∆µ < 30. Negative values of the chemical potential correspond to extremely low ATP concentrations. In particular, the "homeopathic limit" is reached at about ∆µ < −14: At that point there is less than one ATP molecule in an experiment containing one liter of solution.
The reason we still present our results for the entire parameter range −30 ≤ f, ∆µ ≤ 30 is two-fold: firstly, it enables a direct comparison with previous work [11,27]. Secondly, it demonstrates the effectiveness of our algorithm in predicting results that vary over many orders of magnitude. Still, we emphasize that non-trivial results are encountered exactly in the experimentally accessible region, where we consider the model as valid.
A. Velocity, hydrolysis rate, and their quasi tight coupling Fig. 5a reproduces a central result of Ref. [27] concerning the operations modes of kinesin. These modes are defined by the signs of the average currents, i. e. of the velocity c(d) and the hydrolysis rate c(h). However, the resulting phase diagram contains no information regarding their magnitude. Based on the expressions (13) we provide a detailed account on their numerical values in Fig. 5b. Note that these currents vary over about 20 orders of magnitude. This underlines the importance of having analytical expressions to generate the plots. A brute-force numerical approach will either be prohibitively expensive in terms of computer resources, or it will suffer from severe inaccuracies when dealing with this vast range of numerical values.
The analytical expressions for the currents also reveal an interesting relation between the average currents. In Fig. 5d we plot the ratio c(h)/c(d) of the hydrolysis rate and the velocity. Again, access to the analytical expressions for the currents is crucial to determine the ratio. After all, both its numerator and denominator are of the order 10 −18 in the lower left corner of the parameter space.
The most prominent feature of Fig. 5c is that away from the zero-current lines, the ratio of average hydrolysis rate and velocity takes values very close to ±1. Consequently, on average the completion of a cycle yields one mechanical step and one chemical event. We say that chemical and mechanical currents are quasi-tightly coupled. Experimentally, it was found that kinesin hydrolyzes more or less exactly one ATP molecule for each mechanical step [14]. According to the model considered here, quasi-tight coupling is a generic feature that holds Knowing the absolute values of the currents rather than only their signs also allows us to treat kinesin's thermodynamic cycles in more details. In Ref. [27] this discussion was based on the signs of Hill's (excess) cycle fluxes [7]. With the current ratio we interpret the regions shown in the phase diagram Fig. 5a in terms of dominant cycles -at least away from their boundaries: in the upper left and lower right regions the forward cycle F dominates such that average hydrolysis and velocity are directly proportional, c(h) ∼ c(d). The difference between those regions is the angular direction: counter-clockwise completion leads to a forward movement accompanied by ATP hydrolysis, whereas clockwise completion yields backward stepping under ATP synthesis. In contrast, in the upper right and lower left regions hydrolysis and velocity are anti-proportional, c(h) ∼ − c(d): then a counter-clockwise (backward, hydrolysis) or a clockwise (forward, synthesis) completion of the backward cycle F dominates the average dynamics, respectively. This result, which is based on the values of physiological currents, thus complements and extends the one brought forward in Ref. [27].

B. Efficiency of energy conversion
Under physiological conditions, kinesin uses the chemical energy released by the ATP hydrolysis to perform mechanical work. Energy efficiency is one of the most important questions for molecular machines involved in cellular energy conversion [32][33][34], just as it is for macroscopic machines. A framework for a quantitative analysis is based on the notion of conjugate currents and forces from irreversible thermodynamics [23]. Generally, a complete set of conjugate currents c(ϕ i ) and forces E i yields the average dissipation as the bilinear form where σ i denotes the distinct contributions to the entropy production.
Using equations (13) and (14) we find that We see that in the kinesin model, velocity c(d) and hydrolysis rate c(h) are conjugate to the negative pulling force −f and the chemical potential ∆µ, respectively. For a system with two independent contributions to the entropy production, σ = σ 1 + σ 2 , one may define the efficiency of energy conversion in general terms [33]. To that end note that c(σ) is always positive. This, however, does not imply that both contributions c(σ i ) are positive. Indeed, systems act as energy converters only if one of the contributions, say σ 1 , is negative. Then, a (positive) average power outputẆ out := −c(σ 1 ) is sustained by a (positive) average power inputẆ in = c(σ 2 ). Note that c(σ 2 ) is positive and larger in magnitude than c(σ 1 ), because c(σ) = c(σ 1 ) + c(σ 2 ) ≥ 0 must always hold. Hence, the efficiency of energy conversion is defined as It is always positive and smaller than unity.
In the framework of Stochastic Thermodynamics, this efficiency has been studied under various aspects (cf. e. g. Refs. [32][33][34][35]). In Fig. 6, we give the efficiency of energy conversionη for the kinesin model. The regions A-D correspond to different types of energy conversion where the system either acts as a motor (A,C) or a chemical factory (B,D). Outside of these regions both contributions to the entropy production are positive and no energy conversion takes place.
We note the following prediction: for any fixed value of ∆µ in the physiological range, i. e. for 20 < ∆µ < 30, the value of the force at maximum efficiency is around f ≈ 10.5. This suggest that kinesin might be optimized to encounter (elastic) forces of around 5 pN, independent of the ATP concentration. It will be interesting to explore the implications of this result for the collaborative behavior of multiple kinesin molecules involved in the viscous transport of organelles.

C. Diffusion constant and randomness parameter
So far, we have only investigated average currents, which are also available if the steady-state distribution π is known, cf. Eq. (1). Higher order statistics of fluctuating currents cannot be expressed only by means of the stationary distribution, although perturbation expansion exist [36]. The method presented here provides direct access to the (co-)variance of fluctuating currents via Eq. (10b) -and thus without any knowledge of the stationary distribution.
For motor proteins we are mostly interested in the sec- In Fig. 7a we show the diffusion constant in the sixstate kinesin model. Like the average velocity, its values span a range of about 20 orders of magnitude. Under physiological conditions (f = 0, ∆µ ≈ 25) the diffusion constant is about ten orders of magnitude larger than at equilibrium.
A direct measurement of the parameter dependence of D is difficult. An observable that is more easily accessible in experiments is the so-called randomness parameter (sometimes called Fano factor) [15,17,37,38] It is a dimensionless measure of the temporal irregularity of the mechanical displacement. While r = 0 indicates a deterministic motion without any fluctuations, a value of |r| = 1 amounts to a Poisson motor [37]. In Fig. 7b we plot its inverse, r −1 , which is a smooth function. We see that the six-state model predicts Poissonian behavior |r| ≈ 1 in a large area away from the stalling lines. This is in agreement with recent experimental results and theoretical predictions from an alternative model [17]. Our method to calculate the second scaled cumulant and thus the diffusion constant avoids all of the combinatorial complexity of previous approaches [38,39]. Ref. [40] treats drift velocity and diffusion in Markovian models formulated for a periodic lattice in arbitrary dimensions. In the present work the topology of physical space is independent from the structure of the graph representing the topology the model: If a system like a molecular motor moves in more than one spatial dimension, one simply defines a jump observable d i for each physical dimension i. Up to a factor of two, the scaled covariance tensor c(d i , d j ) then equals the diffusion tensor.

D. Response theory
Eq. (15) states that the average velocity c(d) and hydrolysis rate c(h) are conjugate to the mechanical and chemical driving forces −f and ∆µ. Response theory studies the dependence of averaged currents J = (c(ϕ i )) i to the conjugate fields E = (E i ) i . For B independent conjugate current-field pairs, (c(ϕ i ), So-called fluctuation dissipation relations (FDR) relate the response of average currents to their fluctuation statistics. In particular, the Einstein relation relates the mobility of a particle (or its inverse, the friction coefficient) to its diffusion constant. So called Green-Kubo relations express equilibrium transport coefficients by timecorrelation integrals, Eq. (6). [41] In the present context time-correlation integrals are obtained as second-order scaled cumulants c(ϕ i , ϕ j ), and the famous fluctuation relation for the entropy production ensures the validity of the following equilibrium FDR [21]: With analytical expressions for the average currents c(ϕ i ) one can directly calculate their derivatives R ij . Because the correlation integrals c(ϕ i , ϕ j ) are known, our method enables us to probe the non-equilibrium response properties predicted by models of Stochastic Thermodynamics. As an example, we discuss kinesin's mechanical response using the normalized response coefficient The equilibrium FDR (19) implies that T mech (0, 0) = 1. As the transition matrix depends smoothly on the driving parameters (−f, ∆µ) [11], we expect that there will be a one-dimensional curve in parameter space where T mech (f, ∆µ) = 1. Figure 8 depicts T mech . As expected, we see that T mech (f, ∆µ) = 1 holds along two lines originating from the origin, such that a nonequilibrium FDR holds for these parameter values. Remarkably, one of these lines coincides with the stalling line f = f stall (∆µ), i. e. for parameters where the average velocity vanishes.
Another non-trivial feature of Figure 8 is the region where the normalized mechanical response is negative. Note that since c(d, d) = 2D > 0, T mech < 0 implies that ∂c(d)/∂(−f ) < 0. This phenomenon is known as negative differential mobility [42], or more generally, negative differential response (NDR) [18,43]. The kinesin model predicts negative differential mobility for large enough pulling forces beyond stalling, i. e. in situations where the motor walks backwards. Then, by pulling more on gets liess, i. e.the velocity in pulling direction becomes smaller. Interestingly, this feature might already be visible in the experimental data found in Refs. [16,17]. Although we do not expect to see NDR for arbitrarily high pulling forces in real experiments, explicitly looking for it in the region for small superstalling forces seems worthwhile.

E. Model comparison
Direct access to the non-trivial features of physical currents allows us to compare different models in detail, both qualitatively and quantitatively. We start by quantitatively comparing the results for the six-state model ( Fig. 3) with the simpler four-state model (Fig. 4). Recall that the latter is constructed following the same physical arguments as the former (cf. Appendix A for the details). The results plotted in Figs. 5-8 are all derived from c(d), c(h) and c(d, d) = 2D. In Fig. 9 we plot the relative deviations of these quantities between the six-state and the four-state model. Throughout most of the parameter space, they are only a few percent. This is remarkable, because the observables themselves vary over many orders of magnitude. Note that at the boundaries between different operation modes (Fig. 5a), the first cumulants vanish. Hence, we have a divergence in the relative errors unless this happens exactly at the same parameter values in both models.
For the hydrolysis rate c(h) such a divergence shows up in Fig. 9b around (f, ∆µ) (14,14). In principle, this divergence is present wherever c(h) vanishes in the  [11]. We show the relative errors X4/X6 − 1 of the corresponding quantities X4 and X6 calculated in the four-and six-state models, respectively. Throughout the parameter range considered, they are almost everywhere well below 15%. Note that for the average hydrolysis rate the relative error diverges close to the line c(h) = 0 where the hydrolysis rate in the six-state vanishes. Remarkably, this is not the case for the average velocity, where the stalling lines in both models agree exactly. six-state model. In practice, however, the curves of zero average hydrolysis rate c(h) = 0 agree almost perfectly, such that the region where the divergence has an effect is extremely small. For most parameters it is hidden due to the finite (plotting) resolutions, and thus not visible in Fig. 9b. In contrast, the prediction of the stalling forces f stall (∆µ) agrees exactly between the two models: Fig. 9a does not exhibit any singularities. In Ref. [13] we introduced a coarse-graining procedure which preserves the cycle topology of a model. By construction, the first cumulants of all currents agree between the original and the coarse-grained models. Moreover, the relative error in the diffusion constant is comparable in magnitude to what we see in Fig. 9c. These quantitative results emphasize the value of the cycle perspective introduced in Sec III D: In order to construct thermodynamically consistent models, one should think of the physics of cycles rather than focusing only on individual transitions. Finally we compare the six-state model to a general model for molecular motors presented in Ref. [44], which was studied in detail in Ref. [12]. Unlike the six-state model studied fo far, that model features only two states, which correspond to a strongly and a weakly bound configuration. Multiple transition between these two configurations are possible and represent either an active (ie accompanied by a chemical event) or passive displacement along the microtubule. The cycles of the two models are different both in their topology and their interpretation. In particular, the two-state model of Refs. [12,44] has no reference to the "hand-over-hand" stepping mechanism of the forward cycle of Ref. [11] depicted in Fig. 2b. Moreover, the two-state model was fitted to the experiments of Refs. [14,15], whereas the six-state model used the experimental data from Ref. [16] Due to the simple structure of the two-state model, an analytical parameter-dependent expression of the scaled cumulant generating function was found in Ref. [12]. Consequently, analytical expression for the scaled cumulants are known and can be compared to the results ob- 10. Normalized mechanical response T mech for the sixstate model of Ref. [11] studied in the present work (left), and the model from Ref. [12] (right). For experimentally sensible parameters both models predict the same qualitative behavior: The validity of a non-equilibrium Einstein FDR (green curves) at stalling (black curve) as well as negative differential response (region to the right of blue curve) for superstalling forces. Other features (like the overall structure, magnitude of the response) show the same qualitative behavior between the models. Beware, however, the different scaling of the horizontal axes.
tained for the six-state model of Ref. [11], which we used so far. Due to the different nature of the models, we do not expect their predictions to agree quantitatively. This is especially true for parameter values that are far away from the parameters used in the experiments, i. e. for small (or even negative) chemical potentials ∆µ or negative values of the pulling force. However, for experimentally accessible parameters, it makes sense to look for qualitative agreements in the features of the two different kinesin models. Fig. 10 shows the normalized mechanical response Eq. (20) in both models for sensible chemical potential differences (5 ≤ ∆µ ≤ 30) and positive pulling forces. As expected, the models do not agree quantitatively. In particular, the stalling lines are at different positions. However, they show the same qualitative features: the validity of an Einstein FDR at stalling, and the emergence of negative differential mobility just above stalling. Together with the experimental hints from Ref. [16], we consider this agreement as evidence that negative differential mobility is a generic feature of kinesin -a prediction that should be studied by future experiments.

V. CONCLUSION
In the present work, we gave an explicit procedure for the analytical, i. e. fully parameter-dependent, calculation of the statistics of fluctuating currents in Stochastic Thermodynamics. The algorithm applies to any finite Markov model. We focused on its efficiency in exploring the parameter space of models for the motor protein kinesin, while the mathematical background of the algo-rithm was the subject of an accompanying paper [10]. In the following we summarize conceptual and physical insights.
From a conceptual point of view, we find the following points particularly noteworthy: • Our algorithm is efficient. After obtaining the fundamental chord cumulants, Eqs. (10) provide fully parameter-dependent expressions for averages and time-correlations of arbitrary currents.
• Numerical evaluation of expressions is only necessary at the time of plotting. This prevents floatingpoint inaccuracies even in expressions that vary over many orders of magnitude, cf. Fig. 5b or Fig. 7a.
• Having access to analytical expressions allows further (algebraic) manipulation and thus the study of derived expressions, cf. Fig. 5c or Fig. 7b. Taking derivatives with respect to external parameters is needed to explore response properties, cf. Eq. (20).
From a physical perspective, our method allows the systematic comparison of the predictions made by various kinesin models, cf. Sec. III. In particular, we gave a detailed account on (quasi-)tight coupling, efficiency of energy conversion, diffusion and mechanical response for a well-known kinesin model [11] in Secs. IV A-IV D. Moreover, in Sec. IV E we compared these predictions with other models. Regarding the modeling of molecular motors, we emphasize the following: • Current statistics correspond to experimentally observable quantities, like the average motor velocity or its nonequilibrium diffusion constant. Our systematic approach thus extends and unifies previous approaches for calculating these quantities [38][39][40].
• Thinking of stochastic models in terms of its physical cycles is useful. It allows model reduction, cf. Sec. IV E and Ref. [13].
• Independent models predict two interesting nonequilibrium response properties of kinesin: (i) the validity of a non-equilibrium fluctuationdissipation relation at mechanical stalling, and (ii) negative differential mobility for superstalling forces. Both of these predictions lie in realistic parameter regions and can be tested in future experiments.
Modeling the dynamics of molecular motors as random transitions on a biochemical network of states is only one of many appliations of finite Markovian jump processes. The methods established in the present paper apply to any other dynamically reversible model and are easily extended to systems with multiple transitions between states. Thus, they provide a powerful framework to fully explore the physical predictions of any model described by Stochastic Thermodynamics.
are obtained as first-order rate constants κ i j , which are multiplied by concentration and force-dependent factors.
In accordance with first-order rate kinetics, the chemical factor reads For chemical concentrations which are not too high, the non-dimensional chemical potential difference is given by ∆µ = ln K eq [ATP] [ADP] [P] , where K eq = 4.9 × 10 11 µM is the chemical equilibrium constant for the ATP hydrolysis reaction happening at kinesin's active sites. Like in Ref. [11]  For the six-state model the transitions of the forward cycle F = ζ 1 reflect experimental data. We briefly outline how we use the arguments of Ref. [11] for the parametrization of the four-state model shown in Fig. 4. First note that transitions associated to the edges (1 ↔ 3) and (1 ↔ 4) are also present in the six-state model. We thus use similar parameterizations. Transition (3 → 4) combines ADP release at the leading site with hydrolysis (and immediate P release) at the trailing one. In the sixstate model, the same numerical value of the mechanical Values correspond to the experimental data in Ref. [16] as stated in Ref. [11]. The equilibrium constant of the ATP hydrolysis reaction is Keq = 4.9 × 10 11 µM. The parameter θ = 0.65 enters the mechanical factor of the transition rates.
parameter χ i j = 0.15 is used for both of these transitions. We take this as a motivation for using χ The parameters of the remaining transitions are obtained by symmetry. The exception is the first-order constant κ 3 2 , associated with ATP release from the leading head. Similar to Ref. [11] we adjust it in order to account for the Hill-Schnakenberg condition. On the dissipative cycle D = ζ 2 this constraint amounts toσ 2 = 2∆µ and yields κ 3 2 = κ 3 1 κ 1 3 2 κ 1 4 . At this point, we have determined all the parameters of the four-state model while ensuring the physical and thermodynamic consistency with the original six-state model. Fitting yields κ 4 3 = 2.52 × 10 6 , which proves to be a good choice globally, cf. Fig. 9. All model parameters are summarized in Table I.