Analysis of kinesin mechanochemistry via simulated annealing

The molecular motor protein kinesin plays a key role in fundamental cellular processes such as intracellular transport, mitotic spindle formation, and cytokinesis, with important implications for neurodegenerative and cancer disease pathways. Recently, kinesin has been studied as a paradigm for the tailored design of nano-bio sensor and other nanoscale systems. As it processes along a microtubule within the cell, kinesin undergoes a cycle of chemical state and physical conformation transitions that enable it to take ~100 regular 8.2-nm steps before ending its processive walk. Despite an extensive body of experimental and theoretical work, a unified microscopic model of kinesin mechanochemistry does not yet exist. Here we present a methodology that optimizes a kinetic model for kinesin constructed with a minimum of a priori assumptions about the underlying processive mechanism. Kinetic models are preferred for numerical calculations since information about the kinesin stepping mechanism at all levels, from the atomic to the microscopic scale, is fully contained within the particular states of the cycle: how states transition, and the rate constants associated with each transition. We combine Markov chain calculations and simulated annealing optimization to determine the rate constants that best fit experimental data on kinesin speed and processivity.


INTRODUCTION
Kinesin is a molecular motor protein that processes hand-over-hand on microtubules, at each of its 8 nm steps hydrolysing a single molecule of ATP (1,2). Hydrolysis occurs at one of the two catalytic domains (heads) responsible for procession (3). These interact through the neck linker, a polymer of 10-16 amino acids in length connecting the heads. The walk is initiated by ATP binding to one of the heads, causing the neck linker to dock (4)(5)(6)(7)(8). The resulting strain on the second head biases its diffusional search toward the next binding site on the microtubule (9)(10)(11).
Each step in kinesin procession can be modeled as a cycle composed of stages in which kinesin assumes different chemical states and physical conformations (12)(13)(14). This cycle is the starting point in developing a kinetic model for kinesin, in which all information about the kinesin step is contained in the stages of the cycle, in how they transition from one another, and in the rate constants associated with each transition. For this ability to embody information on the mechanochemical gating of the kinesin step and yet being relatively simple, kinetic models are the usual choice for performing numerical simulations of the kinesin walk on a microtubule (15)(16)(17)(18). However, the small length scales and short time scales prevent experiments from resolving each stage of the kinesin step. Therefore the catalytic cycles and resulting kinetic models, must be deduced from what is known about ATP hydrolysis, interaction of kinesin heads with the hydrolysis products, and neck linker physical characteristics.
Our numerical investigation of the kinesin walk aims to reproduce the accepted experimental results for kinesin velocity and run length as a function of ATP concentration and at different experimental conditions with respect to the concentration of ADP and inorganic phosphate (P i ).
The goal of performing an optimization routine involving simulated annealing to fit the experimental data for kinesin speed and processivity is two-fold: firstly, we aim to model the kinesin cycle in a unified way, including all states allowed by experiments. This way we can identify the states that are visited often during the walk (common states) and those which are visited rarely, but may have an impact on the overall performance of the protein walk. Secondly, we can find "better" rate constants: at the end of the simulation, some rates will be very close to the values found experimentally, but others will be off. This means that either they are not contributing to the kinesin walk, or they need to be better investigated by experiments and theoretical models.

Mechanochemical-state model
The first step towards finding a kinetic model for kinesin is choosing which states will be present in the model rate constants and transitions between states follow from this choice. We aim to be as general as possible, utilizing only experimental data to constrain our model.
We start by listing all states that could be included in kinesin mechanochemistry, which can be found by combinatorics: there are 4 different chemical states for each of the two catalytic cores, which can be empty (APO state, also known in the literature as the rigor state (19)) or bound to either ATP, or ADP, or ADP+P i ; and two mechanical conformations of the each head with respect to the microtubule (bound or unbound to the MT). Combining all these possibilities and assuming that the heads are distinct, we obtain 80 possible states. If the catalytic cores are indistinguishable in the stochastic simulation, the number of states is reduced to 48. This number can be further reduced by incorporating the following assumptions based on experimental evidence: (i) kinesin fully detaches from the microtubule and ends its processive run in the ADP-ADP state (20); (ii) the neck linker is in the docked configuration on any head containing ATP or ADP+P i , and is undocked otherwise (9,21,22); and (iii) if both heads are bound to the microtubule, the neck linker cannot be in the docked configuration on the front head due to strain (23). In the third assumption, even though (23) also forbids a two-head bound state with the neck linker undocked in both heads, those states were added due to evidence that some of these states may appear in mechanisms that include backstepping and futile hydrolysis (14). We are then left with a 25-state model that embodies all possible chemical bonding states in the system, taking into account the three essential features of the kinesin catalytic core during the cycle: nucleotide bonding, neck-linker configuration, and microtubule binding. This 25-state model includes all states commonly proposed in the literature for the kinesin catalytic cycle (see, for instance Refs. (15,18,24,25)).
In our stochastic chemical kinetics scheme, the set of discrete states to be used in the model is first enumerated as shown in Fig.1. These states correspond to different nucleotide occupations in each catalytic core, distinct neck linker conformations (docked or undocked to each head) and the various binding states with respect to to the microtubule (one-head, two-head or detached). In Fig. 1(a), we have all the allowed two-head-bound (2HB) states, and in Fig. 1(b), we picture the allowed one-head-bound (1HB) states. We identify as chemical transitions all those indicated by the arrows in Fig. 1. In order to distinguish the mechanical transitions, we double the 8 allowed 2HB states offset by two rows and set Fig. 1(a) below Fig. 1(b), in such way that the 16 states are perfectly stacked, as shown in Fig. 2, left. Figure 2 also includes the detached state (0HB), on the top layer. Labels corresponding to the unique states are shown in Fig. 2, right. Once the arrays of chemical transitions are disposed in layers, we can characterize the mechanical transitions as the ones that cause a change in level (attachment/detachment of one head to/from the microtubule). Such transitions are indicated by the arrows in Fig. 3. In summary, Figs. 1-3 show all possible transitions between the 25 allowed states, where transitions along the "x"-direction are chemical (within a level) and transitions along the "y"-direction are mechanical (interlevel).
The transitions to each state, represented by the arrows in Figs. 1-3, are associated to a set of rate constants {k i→j }, related to the probabilities of transitioning from state i to state j. The rate constants we use to initialize the problem are found in Ref. (26) and references therein (see Table 1, column "Initial"). Figure 4 shows the matrix form of the rate constants, where the states i, j correspond to the header column and row, respectively, by the numerical labeling shown in Fig. 2. A matrix of rate transitions ω i→j can be derived from the rate constants matrix of Fig. 4. We use the relation ω i→j = k i→j I i→j {[X]} from (16), where the indicator function I i→j {[X]}, is unity for every transition that does not involve the capture of ATP, ADP or P i by any of the catalytic heads, and is equal to [X], if the chemical reaction results in attachment of X, where X can be ADP, ATP or P i . Combined, the set of 25 discrete states and matrix of rate transitions based on Fig. 2 contain all information necessary for a kinetic simulation of the kinesin stepping mechanism.

Kinetic Monte Carlo
The sequence of states assumed by kinesin is assumed to be Markovian, meaning that the system carries no memory of states previously occupied, and moving to the next state in the process depends only on the current state of the system. For a Poisson process the time to transition from one state to the other has exponentially decaying statistics and the probability density function of leaving a state i for t > 0 is: where a i = j ω i→j . The total time needed to leave a state i is (27): We use the property of exponentially decaying transition probability density functions for the time to escape a state in kinetic Monte Carlo and Markov theory to estimate the mean run time of the protein and therefore describe the kinesin stepping mechanism.
Kinetic Monte Carlo (KMC) is a stochastic algorithm developed to study the dynamics of chemical reactions (28,29). The algorithm consists in drawing two random numbers r 1 , r 2 in [0,1] before each state transition. The first random number gives us the time to transition ∆t i→j from state i to state j according to the rule: with a i as defined previously. The time to transition is therefore drawn from an exponentiallydecaying distribution. The state j to which the system will transition after a time ∆t i→j is the first state for which the following inequality is true: with a (j) i = j l=1 ω i→l . The simulation runs until the detached state (labeled 25 in Fig. 2) is reached. The advantages and disadvantages of kinetic Monte Carlo stem from the very simple and repetitive nature of the method: whilst we can generate a long chain of states and therefore obtain a detailed description of which transitions occur in each run, meaningful observables such as run length and run time can only be extracted by averaging over many runs, making the method computationally expensive to be used in the optimization.

Markov Chain
The Markovian nature assumed for the dynamical process grants us the use of Markov theory to determine the observables of the system (14,28,30,31). The 25-state kinetic model is discrete, but states can transition at any time t. Therefore, our model corresponds to a continuous-time Markov chain (CTMC), where the time spent in each state is associated to an exponentially decaying probability, as given in Eq. 2. To determine the number of times a state is visited before reaching the detached state, the homogeneous CTMC can be reduced to an embedded Markov chain, treated as a discrete-time Markov chain (DTMC). One major characteristic of the 25-state model is that every state can be visited at a time t until reaching the detached state, when the walk stops. In Markov theory, the states that can transition to another are called transient states, and the states that once attained, cannot transition to any other state are known as absorbing states. For the 25-state kinetic model of Fig. 2, the absorbing state is the fully detached ADP-ADP state, and all others are transient states. Therefore, our embedded Markov chain is an absorbing Markov chain.
The matrix elements of an absorbing Markov chain are the transition probabilities to go from state i to state j, defined as p ij = ω i→j / i =l ω i→l , for i = j, and p ii = 0, except for the absorbing state, where p ii = 1. The transition probability matrix is of the form: where Q is the transition probability matrix between transient states, R is the transition probability matrix to go from the transient states to the absorbing states, I is the identity matrix, and 0 is the zero matrix. From the transition probability matrix Q, we can determine the number of times f ij that state j appears if we start the cycle at state i (using Einstein summation rule): Therefore, the matrix elements f ij of the fundamental matrix F are the average number of times that a state j appears in the chain if the process starts at state i: where q (n) ij is the i, j element of the Q n matrix. For all i, j, we can construct the fundamental matrix F using the infinite sum formula: since q ij < 1 ∀ i, j, the expression above is valid. Once the fundamental matrix is known, we can determine the number of times n jk that a transition j → k occurs if the cycle is started at state i, as the product between the number of times that state j appears and the probability of moving from j to k: n jk = f ij p jk (no Einstein summation rule). To find the average run length of kinesin starting the cycle at state i, x (i) , we define a step matrix S whose elements are: s jk = 0 if the transition j → k does not result in a step; s jk = 1 if the transition j → k results in a step forward; and s jk = −1 if a step backward occurs every time kinesin's state changes from j to k. The average run length is: where ∆x is the 8.2 nm step size. We know from Eq. 2 that the average time τ j to leave a state j is equal to a −1 j = 1/ j =k ω j→k . The average run time t (i) is therefore the sum over all states j of the product of the total time to leave state j and the number of times that state j appears on average, if the Markov process starts at i: Thus the run length and run time are calculated by two independent methods: kinetic Monte Carlo and Markov theory. Markov theory is used within the optimization to calculate the observables numerically. Once the optimization step is finished, we use kinetic Monte Carlo to find error bars associated with the observables.
In Markov theory, we calculate the mean velocity v M C = x / t , from the mean run length x , and mean run time t . This velocity is the quantity calculated by (20). It is important to point out that this speed is not the same as v = x/t , that can be determined from kinetic Monte Carlo. However, in all our calculations the values of v M C and v are within the experimental error bars in (32), which justifies the use of Markov theory within the optimization scheme.

Learning and Optimization Algorithm
Within the 25-state kinetic model, the rate constants found elsewhere in the literature (26) are systematically changed by an optimization process until we obtain an agreement between our numerical results and the experimental data shown in (20,32).
In addition to the states that are part of the kinetic model, the rate constants related to the particular mechanical and chemical transitions connecting these states play a vital role in the determination of the observables of the system, as explained above. Consequently, changes in the rate constants will be followed by changes in the speed and processivity of kinesin as a function of ATP and hydrolysis product concentrations. The rate constants provided in the literature via modeling for the 25-state model (c.f. Table 1, column "Initial" (26)) cannot reproduce correctly the experimental data on kinesin speed versus ATP concentrations at different experimental conditions. The present scheme systematically varies the rate constants in order to train the 25-state model to reproduce a subset of existing experimental results. The model is then tested against "out-of-sample" data to assess its robustness when presented with new data.
We define our objective function Ω as the Euclidean distance between the experimental and predicted measurements: where f exp (y) is the experimental data and f num (y) is from the kinetic model in a given simulation. We further define a normalized objective function Ω n : where N is the number of data points, and: .
Here the variable y corresponds to an independent variable used in the experiment, such as ATP, ADP, and/or P i concentrations. The functions f exp (y) and f num (y) and the expression for Ω n apply to any observable of the system, e.g. speed v or processivity x. In practice, it is convenient to define a combined objective function φ n ({Ω n i,j }), which simultaneously optimizes with respect to multiple observables i (e.g. speed and processivity) under j experimental conditions (nucleotide concentration, etc). The normalized function φ n ({Ω n i,j }) is defined as: with n s equal to the number of experimental conditions, and n o is the number of observables. By construction, at t = 0 the values of the objective function φ n ({Ω n i,j }) lie in the interval [0,1]. For simplicity, from now on we refer to the normalized objective function as simply φ.
Once an initial value of φ is determined for the original set of rate constants, we can start an optimization routine to minimize φ as a function of varying rate constants. The optimization utilizes simulated annealing. The system is started at a high temperature T high with the same initial rate constants used to determine φ. We define β as the inverse of temperature, and assume the Boltzmann constant k B = 1. We run the simulation 500 times, and in each step m, β m is varied in constant steps ∆β, which consists in variable temperature steps that become smaller as T m → 0. At each temperature step, a Metropolis search is repeated 1000 times. The Metropolis step consists of selecting one of the rate constants at random, changing this rate constant, again based on a random procedure, and using a probabilistic criterion based on the current temperature T m and objective function φ to determine whether the change should be accepted or not.
The specific algorithm is as follows: 1. Initialize rate constants as given in Table 1, column "Initial".
2. Choose three random numbers: an integer r 1 ∈ [1,n], where n is the number of rate constants, a real r 2 ∈ [0,1] and an integer r 3 ∈ {1,2}. The first number r 1 is the index of the rate constant we are modifying, k r 1 . The change is done by: so that the new rate k r 1 ∈ [0.75k r 1 , 1.25k r 1 ].
3. All the observables are recalculated for the new set of rate constants {k i } using the Markov chain method.

4.
A new objective function φ new is determined and the acceptance probability is calculated as follows: If the change is accepted, k r 1 → k r 1 and φ → φ new ; otherwise the original rate constants remain the same. Return to 2 until final Metropolis step is concluded.

5.
We record the lowest value of φ given in the 1000 Metropolis steps and the rate constants that generated this value, and use them as a starting point at the next temperature T m+1 < T m . In some cases, we may want to fine tune the fitting, which is done by restarting the simulated annealing at T high but with the values of {k i } that were previously optimized. This prevents the parameters from being trapped at a local minimum and may result in a even smaller value of φ at the end of the second optimization.
6. Once the optimization is finished and the value of φ is low (of the order of 10 −3 or lower), we use the new set of rate constants to determine the observables now through the kinetic Monte Carlo method. For each data point, we run KMC 30 times for different random seeds, and all results are averaged over the realizations.
The optimization is done for 100 different initial conditions, and each one will reach a different solution in parameter space. The solution that gives the best fit is the one we choose to model our system. We discuss our results and observations in the next sessions.

Model Implementation
In order to assess the ability of our modeling to refine the values of the constituents rate constants in an unbiased way, we utilize a train/test machine learning approach. In each numerically optimized experiment described below, we perform the Markov chain/simulated annealing fitting procedure by initializing the entire set of rate constants, performing 100 realizations, and extracting the final "best fit" rate constants within our mechanochemical model that best reproduce a specified set of training experimental data, measured under various [ATP], [ADP], and [P i ] concentrations. Once the optimized rate constants have been determined, they are held fixed, and we compute model predictions for additional experimental data not used in the fit in order to assess, without bias, the quality of the model and thus the computed rate constants.
The experimental data used include velocity vs.
It is important to point out that Ref. (32) shows data for wild-type kinesin and Ref. (20) studies the RK430G kinesin mutant. Kinesin RK430G presents the same speed as wild-type kinesin but its processivity is shorter (20). For our model, the only rate affected for shorter processivity at the same speed is the detachment rate, (k 32 in Table 1). Therefore we can safely assume that all the other rate constants should be the same for both kinesin types. Once we have the optimized set of rate constants (shown in Table 1, column "Optimized"), we run a KMC simulation to obtain the position of the molecule as a function of time at different ATP concentration. The x(t) curves are fitted into a first-order polynomial in t by a least-squares fit, for each independent KMC run (function polyfit in MATLAB). The speed of kinesin for that particular run is simply the slope of the fit. Error bars are obtained from KMC calculations, calculating the standard error of the mean over 30 realizations. It is important to emphasize that our objective function only comprises experimental data on processivity and speed, not constraining the number of ATP hydrolyzed per step. During the kinetic Monte Carlo simulation, after the optimization, we can easily count how many ATP molecules are hydrolyzed every time kinesin steps forward.

Numerical Experiments
Here we explain each experiment done throughout the optimization procedure, starting from curves obtained using the initial assumed rates. Table 2 shows the concentrations of ADP and P i in each numerical experiment performed (training phase).

No optimization:
In this experiment, we use the rate constants compiled from the literature and calculated using the Wormlike chain model for the neck linker (26). The rates are listed in Table 1, column "Initial". We try to reproduce the curves of speed for two cases, both at [P i ] = 0.0 mM and: [ADP] = 0.0 mM (Fig. 5A), and [ADP] = 1.0 mM (Fig. 5B). As observed, the initial rates fail to reproduce the expected curves for kinesin speed.  (20,32). Once the optimization is completed, we obtain a set of optimized rate constants. We then use these new rates to restart the system at a higher temperature and redo the annealing procedure in order to find a lower value of the objective function. As there was no inorganic phosphate in solution during this training phase, the optimized rates cannot correctly predict the speed and processivity data for the experiments with inorganic phosphate in solution. We use this set of rates as the initial set to the following numerical experiment. For the experimental set containing P i in solution, we assume that P i has no effect on the processivity (32) and use the same processivity value for the experiment with [ADP] = 1.0 mM and [P i ] = 0.0 mM. The presence of inorganic phosphate is necessary to constrain the rates related to P i capture, and therefore, will in principle be able to reproduce experiments with nonzero [P i ] during the testing phase. The results of the training phase are shown in Fig. 6, green data (upper right corner).

Testing Phase
After finding the optimized rates and how well they reproduce the data included in the objective function, one should test whether the rates are able to reproduce data that was not previously included in the optimization routine. The test phase is done for the data in (32)  One important result of the optimization is that after Experiment B, the best results in both the training and testing phases recover the ∼ 1 ATP hydrolysis per step forward, as observed experimentally (1,2). Table 2 shows the ratio of ATP hydrolysis per step forward for each of the experiments and the testing phase after Experiment B.

DISCUSSION
Mechanochemical cycles for kinesin have been widely proposed since the first experimental results that allowed a more detailed investigation of particular states of the kinesin walk. Our mechanochemical model, a general network of states, builds upon previously prescribed states and pathways for kinesin walk, based on both theory and experiments.
An early review by Vale and Millikan (33) compares the kinetic cycles of kinesin and myosin by presenting a simplified mechanochemical model very similar to the cycle proposed by Rice et al. (9). The latter was the first work to identify the conformational changes in the neck linker as the force-generating action responsible for the mechanical step. A 1HB start state of empty-ADP was also present in this model, however it was not yet identified as the ATP-waiting state. As experiments extracted more information and gave stronger evidence for a hand-over-hand model of kinesin walk, proposed models became more complex, adding more intermediate states between steps (10, 24) and possible branches (24). Experimental evidence also settled on the existence of certain states of the model for conventional kinesin, such as the ATP-waiting state (34) as the 1HB empty-ADP state, and the ADP-ADP detachment state (20). On the other hand, some states were proven to be forbidden in conventional kinesin due to neck linker strain between the two heads (23). Some theoretical constructs were based largely on kinetic cycles assuming thermodynamic equilibrium fitted to experimental data. Such models were advanced by (15,18,35). The works of Liepelt et al. (17,18,36) pioneered in presenting a clear distinction between which transitions are considered chemical or mechanical, yet treating them at the same theoretical level when building the mechanochemical model. They also introduced a network theory that offers multiple possible pathways for kinesin stepping. The presence of these shortcuts allows for backstepping and futile hydrolysis, such as found by Clancy et al. (25), in which they propose an universal model with branches for a kinesin mutant. Most of the kinetic cycles discussed above are subsets of our proposed network model, as shown in Fig. 7 for the pathways of Refs. (15,18,24,25) (some of these cycles assume states that we have excluded from the model based on previous works (23)). Indeed, further analysis of our results shows that the 25-state mechanochemical model proposed by our work is greatly reduced to a few most-visited states. However there are sporadic visits to rare states, which may become important as they are more or less visited with varying kinesin type and experimental conditions. Thus our model provides a systematic way to analyze the differences in kinesin walk under diverse conditions.
Our results are analyzed once the set of optimized rates is determined via the machine learning algorithm presented in the Results session. The data analysis is done by recalculating the observables (speed and processivity) via kinetic Monte Carlo for 30 different random realizations. Error bars are estimated from the standard error of the mean. A small number of samples is chosen in order to get error bars that are of the same order of magnitude as the ones in Schief et al. (32). Other information that is available through kinetic Monte Carlo are most visited states, time spent in each state and pathways.
States that are frequently visited in our model mostly coincide with states proposed in other kinetic models (12,13,15,24,37). In particular, we find that kinesin not always chooses a specific pathway, but instead takes alternative cycles at each step. Our dominant pathways are shown in Fig. 8. In the absence of hydrolysis products, we found no difference in pathways chosen for different ATP concentrations. In this scenario, about 55% of the time that kinesin takes a step, it does so through pathway 23-11-1-3-5-19-23 in Fig. 8. In 40% of the time, the pathway chosen is 23-11-1-2-4-6-23. One of the main results from Ref. (23) is that kinesin spends about 93% of the time with one head bound to the microtubule, in the ATP-waiting state, at low ATP concentration. We tested our model to check if we observe a similar trend of time spent at the ATP-waiting state. Our results are reported in Fig. 9. Indeed, at low ATP concentration, our results indicate that kinesin spends nearly 90% of the time of the processive run at the ATP-waiting state (93% at 2 µM [ATP]). As expected, as ATP concentration is increased, the time spent at the waiting state shortens. In terms of the actual time spent at the waiting state, our value is about 4 times shorter (≈100ms) than the value found by (23). This is explained by our value of the ATP binding rate being almost 4 times higher than the one reported in (23). Both rates are within the accepted rates of 1-6 µM −1 s −1 (38)(39)(40)(41)(42), and the discrepancy can be explained by the experimental conditions and systems (single molecule vs. gliding motility assay) studied. In fact, most of our rates are well within accepted values, even though they are allowed to vary without constraints. On the other hand, this also results in a few of our rates falling outside the range of experimental values. One advantage of this method is that we are able to determine rates which are not able to be found experimentally, especially those associated to very fast events such as kinesin binding/unbinding to the microtubule.
The change of speed and processivity of kinesin with added ADP is understood by analyzing the most common pathway that the protein undergoes during procession. Without the addition of extra ADP to the experiment, the waiting state of kinesin can only bind ATP, and the cycle of hydrolysis-stepping takes place continuously. As more ADP is added to the assay, the waiting state can bind either ADP or ATP. Whenever it binds ADP, it goes into an ADP-ADP state with one head bound to the microtubule, which is the state that leads to detachment. As this state is more visited, over time the chances of detaching becomes higher thus reducing processivity. Similarly, every time ADP binds to the waiting state instead of ATP, the protein must wait for ADP to unbind before ATP is free to bind and hydrolyse, thus slowing down the walk and reducing the speed.
The simulated annealing method has proven to fit well experimental data by finding new rate constants for the model chosen. The fact that our model is able to recover the rate of 1 ATP per step forward without constraining to it is a major strength of the procedure. We are also able to run very fast simulations by choosing to use Markov chain calculations within the optimization scheme, which are much simpler to solve in each Metropolis step as opposed to running kinetic Monte Carlo (Markov chain requires only a matrix inversion whilst KMC would require hundreds of independent realizations to calculate the average values of the observables). Furthermore, our finite Markov model does not require solving of Master equations to compute observables. A steady-state solution is also not required as kinesin is not allowed to rebind to the microtubule after it reaches the detached state, thus the verification of detailed balance is only necessary if we assume that the system operates near equilibrium. As we chose to solve our model with a minimum of assumptions, we do not constrain our rates to obey equilibrium conditions and therefore most of our transitions are practically irreversible, with reverse rates to some reactions close to zero, characterizing a non-equlibrium system. A difficulty that arises from running stochastic simulations for optimizing rate constants is that is it not always easy to choose which set of rates takes us to the best local minimum and still represent a physical system. Thus one must be careful to analyze the results by comparing the new rates to previously determined experimental rates, and constraining the simulation to more experimental data. Moreover, to provide a more accurate quantitative description of rates, further investigation may be required, such as a sensitivity analysis of the optimized rate constants. Nonetheless, our method is successful in providing a qualitative picture of changes in pathways as experimental conditions change, as we have shown in this work, and enables systematic studies of different kinesin types.

CONCLUSION
In summary, we have successfully implemented a novel simulated annealing method to optimize kinesin catalytic cycle rate constants by constraining kinetic Monte Carlo and Markov chain discrete state simulations to yield average processivity and velocities from experiment. The modeling yields a general network for the kinesin catalytic cycle, capable of making quantitative, verifiable predictions for independent experiments not included in the rate constant optimization procedure. The new discrete state model can be further extended to include additional states for modeling mutant kinesins, in order to elucidate the critical role of neck linker docking in kinesin mechanochemistry. The kinetic model also provides a potential link to the atomistic scale, by highlighting the discrete states of the kinetic cycle whose rate constants are maximally sensitive under optimization to environmental, structural, and chemical changes.      (20,32). Green solid lines: Numerical data obtained from initial rate constants for the model depicted in Fig. 2 via Markov chain calculations. Red dashed lines: Numerical data obtained from optimized rate constants for the model depicted in Fig. 2 via Markov chain calculations.