 Research article
 Open Access
An exact approach for studying cargo transport by an ensemble of molecular motors
 Donatello Materassi^{1}Email author,
 Subhrajit Roychowdhury^{2},
 Thomas Hays^{3} and
 Murti Salapaka^{2}
https://doi.org/10.1186/20461682614
© Materassi et al.; licensee BioMed Central Ltd. 2013
 Received: 7 December 2012
 Accepted: 22 October 2013
 Published: 16 November 2013
Abstract
Background
Intracellular transport is crucial for many cellular processes where a large fraction of the cargo is transferred by motorproteins over a network of microtubules. Malfunctions in the transport mechanism underlie a number of medical maladies.
Existing methods for studying how motorproteins coordinate the transfer of a shared cargo over a microtubule are either analytical or are based on MonteCarlo simulations. Approaches that yield analytical results, while providing unique insights into transport mechanism, make simplifying assumptions, where a detailed characterization of important transport modalities is difficult to reach. On the other hand, MonteCarlo based simulations can incorporate detailed characteristics of the transport mechanism; however, the quality of the results depend on the number and quality of simulation runs used in arriving at results. Here, for example, it is difficult to simulate and study rareevents that can trigger abnormalities in transport.
Results
In this article, a semianalytical methodology that determines the probability distribution function of motorprotein behavior in an exact manner is developed. The method utilizes a finitedimensional projection of the underlying infinitedimensional Markov model, which retains the Markov property, and enables the detailed and exact determination of motor configurations, from which meaningful inferences on transport characteristics of the original model can be derived.
Conclusions
Under this novel probabilistic approach new insights about the mechanisms of action of these proteins are found, suggesting hypothesis about their behavior and driving the design and realization of new experiments.
The advantages provided in accuracy and efficiency make it possible to detect rare events in the motor protein dynamics, that could otherwise pass undetected using standard simulation methods. In this respect, the model has allowed to provide a possible explanation for possible mechanisms under which motor proteins could coordinate their motion.
Keywords
 Molecular motors
 Rare event detection
 Markov models
Background
The behavior of motor proteins is relatively well characterized when one motor protein is involved in the transport of a cargo. Indeed, it is possible to monitor the motion of a single molecular motor under highly tunable experimental conditions and obtain measurements with sufficiently accurate spatial and time resolution [1–3]. The resulting experimental data has led to many theoretical descriptions of motorprotein mechanisms which take into account the complex mechanochemical processes involved and yield insights into transitions between the multiple conformational states possible [4].
In vivo, often, an ensemble of molecular motors is responsible for the transport of a common cargo [5, 6]. In vitro and simulation studies where multiple motors are involved in transport have provided unique insights into features of a common cargo being transported by many motors (see for example, [7, 8]).
The dynamics when multiple motors transport cargo can be considerably more involved where a number of significant questions remain open. For example, it is not yet clear when and if motors synchronize their behavior, whether they move independently and whether they are antagonistically engaged in a “tugofwar” [6, 9]. Despite major improvements in instrumentation and techniques, understanding behavior of multiple coupled motors remains extremely challenging. The main difficulty is the substantially higher spatial and temporal resolution needs imposed by the fractional motion of the cargo and the increased number of possible transitions between conformational states [8, 10]; possibilities introduced by the multiplicity of motors carrying a single cargo.
The available detailed characterization of how single motors transport cargo can be leveraged to develop models that describe how multiplemotors coordinate the motion of a common cargo. Indeed, using single molecule experimental data, accurate descriptions on the probability that a motor takes a step and its dependence on environmental factors such as temperature and ATP concentration, are reported in [11–13]. Similar estimates on the attachment and detachment rates of molecular motors to and from a microtubule can be found in [13–15]. A model that describes how multiple motors carry a common cargo can be obtained by using the information on single motor protein behavior and by introducing the coupling of the individual motorproteins via the dynamics of the shared common cargo. Using MonteCarlo simulations on such a model [7], reported novel insights into the behavior of kinesin motors, such as, a smaller velocity of transport of cargo when carried by multiple motors as opposed to a single one, and a dependence of the expected runlength on the stiffness of the motor linkage. While MonteCarlo techniques form an important set of tools, they involve a tradeoff between the accuracy desired and the computational effort needed. As a consequence, important features of the dynamics, especially if associated with rare events, can be missed. This aspect takes particular significance in the study of biological systems, where pathological behaviors are caused or triggered by events which are improbable under normal conditions but occur with significant adverse impact.
Existing approaches have utilized models with simplifying assumptions that can be treated analytically or semianalytically in order to understand the basic features of the coordinated motion of motor proteins. For example, in [16] meanfield theory is applied for analyzing large ensembles of motors, whereas, in [17] the cooperative transport of cargo realized by two motor proteins is studied in order to identify distinct operational regimes. In [14] apart from providing estimates of attachment and detachment rates of motors to microtubules, analytical dependence of runlength on the number of motors involved in the transport of a common cargo is obtained.
In this article, we present a general methodology which determines the probability distribution function of various motor behaviors. This different approach provides several advantages over MonteCarlo simulation based methods. In our method the probabilities of outcomes are determined exactly, unlike MonteCarlo simulation based methods; however, our method does not sacrifice the detailed description of the system possible with MonteCarlo simulations. Our strategy is particularly well suited for characterizing rareevents that take prohibitive number of simulations in a MonteCarlo setting. Moreover, in the new framework, delineation of the detailed causes of an observed functionality is straightforward (which involves a simple step of identifying states that are associated with the observation and analyzing these states). At the same time, our model has a high level of accuracy and detail. Compared with other analytical studies, such as the ones previously reported, a larger number of motors can be studied. In [18] and [17] the study is limited to only two motors and certain simplifying assumptions are often made (i.e. the aggregation of microstates with same energy in [18]). In [18] a stochastic model that takes into account only the number of motors engaged on the microtubule is adopted in order to understand the level of coupling among two motor proteins carrying a common cargo. In [16] groups of more than two motor proteins are studied. Related work [14] alluded to earlier analyzes the runlength, average velocity, steady state distribution of bound motors and effects of load force on velocities. In both, the meanfield approach of [16] and the approach in [14], the proteins are not individually modeled anymore (for example, it is assumed that the load is equally shared on all the engaged motors). Under the methodology described in this paper, each motor is individually modeled and analytical or semianalytical results can still be provided. Thus, more accurate conclusions on how the interaction between multiplemotors affects a transportation modality can be reached.
The article develops a Markov model, where the number of motors at any particular location on the microtubule lattice form states, and such a state determines the location of the common shared cargo. Here the transition probabilities between states can be derived from studies on single motorprotein based transport. The physics of the system is utilized to project the resulting infinite dimensional model onto a finite dimensional one. We show that the finite dimensional model, apart from the benefit of increased computational tractability, has other important features such as the existence of a unique steadystate probability distribution. Furthermore, we demonstrate that the probability distribution of the projected model can be used to answer most of the biologically relevant queries on transport modality. In particular, probabilities of rare events and the related mechanisms can be unraveled. The capabilities of the methodology are tested with existing data and via extensive MonteCarlo simulations. These features can significantly ease the computational burden as well as provide unique insights into transport modalities.
Methods
Here we provide a methodology for analyzing the dynamics of an ensemble of motor proteins carrying a single cargo on a microtubule lattice. Each individual motor behavior is described stochastically: it can detach or attach to the microtubule and take steps on the filament according to prescribed probabilities that are governed by specified transition rates. The derived stochastic model provides an intuitive representation of the physical system, but, being infinitedimensional, is not tractable and provides no general guarantees on the existence of stationary steady behavior. This impasse is overcome by building an alternative, and effective, Markov model with the advantage of being described by a finite number of states. In this model only the information pertinent to the relative configuration of the motorproteins is incorporated where the relative positions of motor proteins with respect to each other determine the state. The evaluation of the probability distribution for all these possible arrangements can be determined by computing the exponential of a matrix with a dimension that is dependent on the number of arrangements. We show that the number of states does not become excessively large and that the solution via matrix exponential is viable, allowing a direct way to compute the probability distribution of motor arrangements. Furthermore we show how quantities of interest such as, average cargo runlength, average number of engaged motors and average speed of the cargo can be derived, from the determined probability distribution on the relative configurations.
We instantiate the methodology to the case where cargoes are transported by multiple kinesin motors. Despite being specific to these molecules, most of these strategies can be extended or adapted to other classes of motor proteins and also to model a cargo transported by multiple species of proteins, as well.
Description of the system and main modeling assumptions
where k_{ el } is the stiffness of the linkage. If the linkage is stretched beyond a certain stalling force F_{ s }, the motor can not take any forward step. We remark that F_{ s } is typically measured in order to quantify the number of motors that are actively pulling a cargo. Backward steps are neglected in the model and the motors are irreversibly bound to the cargo particle. A motor head that is attached to the microtubule has a certain chance of detaching from it, while a motor head that is not attached has a certain chance of binding to the microtubule. An unbound motorprotein can bind to the microtubule at a location only when it is within a distance l_{0} of the cargo. Thus a floating motor binds to the microtubule without stretching its linkage. The cargo is subjected to a constant load F_{ load } that opposes the motor motion. The cargo position is described in probabilistic terms by a Gaussian distribution with variance σ_{ th } and truncated on the interval [−3σ_{ th },3σ_{ th }]. The mean position of the cargo x_{ eq } is the equilibrium position determined by the load F_{ load } on the cargo and forces exerted by the motors through their linkages. The effect of thermal fluctuations is incorporated into the probabilities of cargo position by determining the variance parameter σ_{ th } of the cargo position in a steady state situation.
When a motor steps forward or detaches, the probability distribution of the position of the cargo is assumed to reach a new distribution with negligible transient. Thus we assume that the time scale of the cargo dynamics is much faster than the rate at which motor configurations change. The system is assumed to be spatially invariant: its stochastic behavior does not change if the motor ensemble and the cargo shift to a new position along the microtubule. Finally, if, at any time, there are no motors engaged with the microtubule, the cargo is assumed to be “lost” which forms the stopping criterion for the stochastic model.
In the model, it is assumed that multiple proteins could share the same location on the microtubule, even though the motor proteins actually bind to physically different areas of the cargo macromolecule. The motivation and justification for this assumption are provided later.
We also observe that the spatial invariance hypothesis translates into an immediate condition on the transition rates, namely that ${\lambda}_{\mathit{\text{abs}}}({Z}^{\prime},Z)={\lambda}_{\mathit{\text{abs}}}({\rho}^{\alpha}{Z}^{\prime},{\rho}^{\alpha}Z)$ for any integer α. This condition, along with the presence of a stalling force for the motors, is used to arrive at an effective finitedimensional Markov model.
Derivation of an effective finitedimensional Markov model
The representation of an ensemble of motors as a biinfinite sequence allows one to describe the system in a rather intuitive manner and highlights the similarities with a Gillespie model for the purpose of stochastic simulations [19, 20]. However, such a model is illsuited for an exact analysis because of its infinite dimension. A finite dimensional model can be obtained by aggregating (or projecting) states of the infinite dimensional model into “macrostates”. In general, this approach leads to the loss of the Markov property. However, in the following we provide a projection of the infinite states of the original model on a finite set in such a way that the Markov property is preserved. This allows us to pursue an exact analysis and determine explicit formulas for the computation of biologically relevant quantities.
This intuitive string representation provides the relative configuration which characterizes how various motors carrying the cargo are positioned with respect to each other.
We make the following observations:

Strings representing relative configurations can have arbitrary length.

Two different absolute configurations of the motors, ${Z}^{\prime}$ and Z, on the microtubule may have the same relative configuration if Z is a “shifted version” of ${Z}^{\prime}.$ Two absolute configurations have the same relative configuration, if and only if the relative distances among the engaged motors of the ensemble are the same. This defines a class of equivalence on absolute configurations: two absolute configurations belong to the same equivalence class if both have the same relative configuration.

From a relative configuration we can obtain the relative positions of the motors, but not their absolute positions on the microtubule lattice.
 1.
An ensemble contains $\overline{m}$ molecular motors (which is the number of motors attached to the cargo)
 2.
Motor linkages are elastic springs with constant k _{ el } and rest length l _{0}
 3.
There is constant load F _{ load } on the cargo
 4.
The stalling force is F _{ s }
 5.
An unattached motor can attach to the microtubule only to locations that are within distance l _{0} from the cargo center of mass (the attachment occurs a locations that are close enough not to stretch the linkage)
 6.
All motors are attached at the same location on the cargo and multiple motors can share the same microtubule location.
The last assumption is introduced for the following reason. From a mathematical perspective, there is no loss of generality on assuming that all molecular motors are bound to the same cargo location. Indeed it is possible to apply a coordinate change to each motor’s position whereby all motors are attached at the same location on the cargo. With this assumption we have to allow for multiple motors to be attached to the same microtubule location, as, identically stretched motors that are physically attached to the cargo at different locations get mapped, in the new coordinate system, as being attached at the same location on the cargo and the microtubule.
where ⌈·⌉ represents the ceiling function. The main intuition on how the various factors in (3) contribute follows from the stall condition on the motors, where, a motor cannot step forward if it experiences a force greater than the stall force F_{ s }. For example, $\frac{\overline{m}{F}_{s}{F}_{\mathit{\text{load}}}}{{k}_{\mathit{\text{el}}}{d}_{s}}+1$ is the maximum distance between the rearguard and vanguard motor possible, beyond which motors stall, $\frac{6{\sigma}_{\mathit{\text{th}}}}{{d}_{s}}$ accounts for the thermal noise contribution, whereas, $\frac{2{\ell}_{0}}{{d}_{s}}$ accounts for the possibility that motors are within a distance 2ℓ_{0} where the motors are not stretched at all.
Let ${s}^{\left(\mathit{\text{max}}\right)}:=\frac{\stackrel{\u0304}{m}{F}_{s}{F}_{\mathit{\text{load}}}}{{k}_{\mathit{\text{el}}}}+2{\ell}_{0}+{d}_{s}$. It follows that if none of the motors are stalled then the extent S≤s^{(m a x)}−d_{ s }.
Now we can assert that if the extent was less than or equal to s^{(m a x)} then for any subsequent change in the configuration, the extent will still remain less than s^{(m a x)}. Indeed, consider the case where the current configuration is such that the extent S≤s^{(m a x)}. There are two possibilities for the current configuration (a) the vanguard motor is stalled in which case the extent can only decrease in any subsequent change in the configuration as the vanguard motor cannot step forward and the rearguard motor cannot step backwards (b) the vanguard motor in the current configuration is not under stall in which case the extent S≤s^{(m a x)}−d_{ s }. In any subsequent change the only means to increase the extent is when the vanguard motor takes a step with a stepsize d_{ s } where the extent still remains bounded by s^{(m a x)}. Thus we have shown that if the extent of an absolute configuration is smaller than a bound s^{(m a x)} then for all future configurations this bound is respected.
Notice that physically different simple events can give rise to the same transition in the symbolic dynamics of the strings. For example, from the string MM it is possible to reach the string M because of a detachment of either the vanguard or the rearguard motor.What has been achieved so far is a projection of the model dynamics from the set of absolute configurations (Z) with infinitely many elements to a space of relative configurations (σ) with finitely many configurations. We denote the projector operator as σ=Π^{(e)}(Z) where the absolute configuration Z has a relative configuration σ.
In general projections do not preserve the Markov property of a model. However, in this case, we can show that the dynamics on the string space still maintains the Markov property. More importantly, the transition rate λ_{ rel }(σ^{′},σ) from one string σ to another string σ^{′} can be meaningfully defined and can be computed from the knowledge of the rates ${\lambda}_{\mathit{\text{abs}}}({Z}^{\prime},Z)$ of the original Gillespie model. We now determine ${\lambda}_{\mathit{\text{rel}}}({\sigma}^{\prime},\sigma )$.
where the first three equalities have been explained before, the fourth follows from Bayes’rule and the fifth is evident. The sixth equality uses translation invariance where the absolute configurations at t and t+Δ t are both shifted by ρ^{−α}, the seventh follows from the fact that the set $\{{\rho}^{(\beta \alpha )}:\beta \in \mathcal{I}\}=\{{\rho}^{\beta}:\beta \in \mathcal{I}\}$ where α is fixed and β is any integer (with denoting the set of integers).
where the $min$ operator has been introduced just to obtain a term that formally depends on σ and ${\sigma}^{\prime}$ only.
where represents the set of all the possible N relative configurations of motorproteins. Thus, the Master equation does hold in terms of the transition probabilities and this implies that the underlying model that governs the dynamics of relative configurations is indeed Markov.
By enumerating the strings σ_{1},...,σ_{ N } that represent relative configurations, we let P_{1}(t),...,P_{ N }(t) represent the probabilities of having the system in each one of the string configurations and define, the probability vector P(t)=(P_{1}(t),...,P_{ N }(t))^{ T }. Using the expressions of the transition rates λ_{ rel }(σ_{ j },σ_{ i }) and Equation (7) it can be shown that
where $exp\left(\mathcal{At}\right)$ is the matrix exponential.
In the specific of kinesin motors, for realistic values of the system parameters and number of motors ($\overline{m}\le 8$), the dimension of A is in the order of 10^{5}−10^{7}, making the problem of computing $exp\left(\mathcal{A}\right)$ manageable for a standard desktop computer. For more complex scenarios (i.e. multiple species of motor proteins or larger ensembles) the problem is still tractable using computer clusters or supercomputers
Determination of biologically relevant quantities
In the previous section, starting from an infinite dimensional model that describes the system dynamics, we have defined a finite dimensional model that keeps track of the relative distances among the motors of the ensemble. The effectiveness of this finite dimensional model is given by the fact that biologically relevant quantities of the system can be computed using explicit formulas without taking recourse to Monte Carlo simulations. Indeed, the probability distribution P(t) of the different configurations provides detailed information about the system, since it provides the probability associated with every specific relative arrangement of motors on the microtubule. Once the probability of having a certain pattern of motors with all the associated relative distances is known, it is possible to determine many quantities of biological interest for the system. In the following, we provide the expressions of certain biologically relevant quantities, as obtained from our finite dimensional model. They will be considered for the validation of the methodology and in the discussion of novel results.
Average number of engaged motors
where M(σ) represents the number of symbols ’M’ in the string σ.
Average velocity and average runlength
We first note that the change in the equilibrium position of the cargo is translation invariant. That is if the initial and the final absolute configurations are translated by the same amount then the change in the cargo position remains unaltered. Thus $d({Z}^{\prime},Z)=d({\rho}^{\alpha}{Z}^{\prime},{\rho}^{\alpha}Z)$ for any absolute configurations Z and ${Z}^{\prime}$.
in order to obtain a right hand side that is formally a function of σ and ${\sigma}^{\prime}$ only.
An important quantity that can be experimentally measured in experiments is the expected runlength of the motors, that is the average length traveled by the cargo/motor complex before movement is arrested or the motor detaches from the microtubule lattice. The average runlength can be determined from the knowledge of the probability vector P(t) on relative configurations.
Distribution of step length
where χ={x:μ^{(x,t)}(x,t)≠0}. This formula provides an exact computation of the distribution of the step size from the model parameters without relying on histograms obtained from Monte Carlo simulations.
Results and discussions
Methodology validation
While the methods developed in this article can be applied to study ensembles of motor proteins of any class, validation will be presented on kinesin motors. We first derive the transition rates λ_{ abs } between absolute configurations for an ensemble of kinesin motors.
Obtaining transition rates on the absolute configuration space
The determination of transition rates are based on experimental data and theoretical considerations, where, rates from singlemotor experiments are used to derive transition rates in the case where an ensemble of motors are involved in transport. Most of the modeling assumption are the same as made in [7] with minor differences which are described next.
Probability of stepping under a force Ffor kinesin
Following [11], MichaelisMenten dynamics predicts a ATP hydrolysis rate equal to k_{ cat }[A T P]/([A T P]+k_{ m }), where k_{ m }=(k_{ cat }+k_{ off })/k_{ on }. In addition, the free head of the motor is assumed to bind to the microtubule location with a defined probability (or efficiency) ε. In this scenario, the probability P_{ step } of stepping
where the term z_{ k } represents the number of motors in the kth location (the transition rate is proportional to the number of motors in the location) and the term a_{ k } is the position of the kth location.
Probability of detachment
Probability of attachment
Experimentally, it is found in [21, 22] that the probability of a kinesin motor attaching to the microtubule is P_{ att }≃5s^{−1}. If the motor is linked to the cargo, it is assumed that it attaches to the microtubule without stretching its linkage. Thus, the only admissible locations of attachment are the locations at a distance from the cargo that is less than l_{0}. They are also assumed all equally likely.
Numerical parameters
The numerical parameters that we have considered in our analysis of KinesinI ensembles, when not otherwise specified, are k_{ cat }=105 s^{−1}, k_{ on }=2·10^{6}M^{−1}s^{−1}, k_{0o f f}=55 s^{−1}, [A T P]=2·10^{−3}M, F_{ s }=0.006 n N, d_{ s }=8 n m, d_{ l }=1.6 n m, δ_{ l }=1.3 n m, A=107, B=0.029 μ M, T=300K k_{ el }=0.32·10^{−3}n N/n m. All these parameters are the same used in [7] and have been experimentally determined.
Using these parameters an upperbound on s^{(m a x)} (see Equation (3)) on the extent of any relative configuration is found to be 320 n m, for ensemble of at most 4 motors. This extent is rather large given that the length of a Kinesin molecule is in the hundreds of nanometer range. We remark that the s^{(m a x)} is an upper bound on the possible extent. Thus there are avenues to be explored where a smaller extent can be assumed. We enumerated the finite number of relative configurations as σ_{1},…,σ_{ N } and determined transition rates λ_{ r e l }(σ_{1},σ_{2}), from a relative configuration σ_{1} to another relative configuration σ_{2}, using Equation (6).
As in [7], we have considered ensembles of at most four motors, and, we computed the probability vector P(t) exactly. We remark that P(t) depends on the initial probability P(t_{0}), as shown in Equation (9). In all our computations for validation purposes we have assumed the same initial probability distribution P(t_{0}) that is used in [7].
Validation of the average velocity and average runlength
Average Runlength: In [7], in one scenario (Model A) the authors neglect the effect of thermal noise, and in another scenario (Model B) they introduce a dynamic model for the Brownian motion of the cargo. We have performed the same kind of separated analysis following our approach. First, we have fixed the variance of the cargo position σ_{ th } to zero, making our model analogous to “Model A” in [7].
For the initial distribution we consider that at time t=−1 s e c exactly one motor is attached to the microtubule and that the cargo is not being lost before time t=0. In the time interval [−1 s e c,0 s e c] the motors behave as usual. The probability distribution of their configurations at time t=0 is the initial probability distribution for all our simulations. This initialization is similar to the one described in [7].
We find a practically exact quantitative agreement of our exact results and the one based on Monte Carlo simulations as presented in [7] which correspond to a coarse grid resolution of 1 pN, when there is no noise (compare Figure 5(a) with Figure 5(a) in [7]). In particular, for all possible sizes of the ensembles, we find runlength curves that are monotonically decreasing with higher loads. In a similar manner a near quantitative agreement is found when noise is present (compare Figure 5 with Figure 5(b) in [7]).
Discussion
The methodology developed in this article allows one to determine how the probabilities of different relative arrangements of molecular motors on a microtubule change over time.
This information is contained in the probability vector P(t) (see Equation 9). For physical values of the model, the number N of arrangements is limited and allows its direct computation. The knowledge of P(t) offers, from a biological perspective, detailed information about the system. In fact, except for the absolute position of the motor ensemble on the microtubule (that is lost in the string representation), the information about the system is completely preserved. When P(t) is known, it is possible to determine, via explicit formulas, many quantities of biological interest, such as the average runlength of the ensemble, the average number of engaged or active motors, the average instantaneous velocity at which the cargo moves and the probability distribution of the step sizes observed in the cargo motion. Contrary to other methods, the final accuracy of the results does not depend on any specific simulation technique or on the number of stochastic simulations that are performed. Also, the method is extremely efficient: even for practically sized ensembles with $(\overline{m}\le 8)$, results can be computed on a standard desktop computer and general purpose software. For a large range of physically meaningful values of the parameters, the number N of possible string configurations is in the order of about 10^{5}−10^{7}. Furthermore, the matrix , that defines the dynamics of the ordinary differential equation to be solved, is sparsely populated. The manageable dimension of the system state and the high sparsity of make the computation of the exponential of feasible even with a limited amount of memory. As evidence for this, all the results shown in this article have been obtained using a machine equipped with a quadcore processor and 4 G b memory (the algorithm was implemented using MATLAB, TheMathWorks, Natick, MA). For the study of larger groups of motor proteins, the adoption of computer clusters or supercomputers still remains a viable solution.
Presence of a steady state
The probability vector P(t) for the different arrangements of motors of the ensemble depends on the initial probability P(t_{0}), as shown in Equation (9). A fundamental question is whether starting from any arbitrary initial condition, after a transient period, the motor ensemble eventually behaves according to a fixed probability distribution which does not depend on the initial distribution of motors. A property of this kind would justify the experimentally observed robustness of the system. Furthermore, it would make it possible to determine the generality of certain observations, independent from the initial distribution P(t_{0}).
The probability of having no engaged motors on the microtubule slowly converges to 1 as t goes to infinity. This corresponds to an intuitive fact: the loss of the cargo is an irreversible event for the system and, sooner or later, it is to be expected that all motors will detach from the microtubule. According to the model formulation, the loss of the cargo is always the final event and, as such, it trivially represents the only steady state condition of the system. However a nontrivial, and biologically more meaningful, notion of “steady state” can be introduced. Figure 7(b) shows the conditional probability distribution of the number of motors engaged on the microtubule at time t, given that at least one motor is engaged. This conditional probability distribution converges to a nontrivial distribution. Thus, under the assumption that the cargo has not been lost, the number of motors reaches a probability distribution that does not depend on its initial condition. This holds not only for the probability distribution of the number of engaged motors, but, more generally for the probability distribution of relative configurations, as we will show at the end of this section.
In other words, under the hypothesis that the cargo has not been lost, the relative arrangements of the motors on the microtubule reach a stable (conditional) probability distribution $\Pi \in {\Re}^{N1}$. The determination of Π, once P(t) is known, is quite straightforward and can be obtained directly from the definition of conditional probability. The knowledge of this steady state Π provides key insights into the behavior of a group of motors.
Data in Figure 8(ab) provide the following insight. The rearguard motor is always assumed as a reference and the xy axes represent the relative distances of the first and second motors of the ensemble from the rearguard one. Thus, each point xy represents an arrangement of the ensemble. On the z axis we report the probability of that specific arrangement. The probabilities of configurations with less than three motors engaged are not reported in the two figures because that would make the visualization difficult. What is important to notice is how the presence of either a low load or a high load leads to two different steady state situations. Thus under a low load the motors tend to spread out almost uniformly, instead, under a high load, a certain pattern of configurations emerge as being more likely. The most likely configurations lie along the diagonal x=y, with a prominent peak around the origin x=y=0. This means that, under a high load, eventually it is more likely to find all the three motors clustered together (represented by the peak at the origin). The high frequency of configurations along the x=y diagonal suggests that it also likely to find two close leading motors with the third one lagging behind. Observations like this would be difficult to obtain using Monte Carlo simulations. Instead, an exact computation of the probabilities allows to infer these characteristics of the motion in a comprehensive manner.
Existence and uniqueness of the conditional stationary distribution
where

each column of the matrix sums up to zero

the upper triangular structure of the matrix derives from the particular way we have reordered the states in accessible from σ_{1} and not accessible from σ_{1}

the top left entry A_{1,1}+A_{N,1} derives from having removed the state σ_{ N }=∅ since we are looking for the conditional distribution of the relative configurations given that the cargo has not been lost.
The bottom right block A^{(22)} is such that (A^{(22)})^{N−1} is a strictly diagonally dominant matrix, since it is possible to reach σ_{1} from each of the states ${\sigma}_{{N}_{a}+1},\dots ,{\sigma}_{N1}$. This implies necessarily that ${\Pi}_{{N}_{a}+1}={\Pi}_{{N}_{a}+2}=\dots ={\Pi}_{N1}=0$. Instead, the top left block of the matrix is an irreducible matrix (in the associated graph each state can reach any other one passing through σ_{1}) implying that the elements ${\Pi}_{1},{\Pi}_{2},\dots ,{\Pi}_{{N}_{a}}$ and are uniquely determined and strictly positive. Thus, there is a unique stationary conditional distribution given that at least one protein is attached to the microtubule.
Enabling finer analysis
As our method yields an exact probability distribution, it facilitates a finer analysis. For example, the dependence of the average run length on the load, under the presence and absence of thermal noise, is of interest. In [7], Monte Carlo methods yield the dependence, where the average runlength is obtained with respect to load in steps of 1 p N. With our method it is straightforward to obtain the exact values at this force resolution. However, we can obtain this dependence at a finer force resolution of 0.2 p N. Using the finer scale, we noticed peaks in the runlength curves for $\overline{m}>2$, that are not evident at the coarse resolution of 1 p N. These peaks correspond to loads that are multiples of the stalling force F_{ s }. These peaks, ascertained by our method provides the following insight. Let us consider the curve corresponding to $\overline{m}=2$ for simplicity. When only one motor is engaged and F_{ load } is close to (but less than) the stalling force F_{ s }, the probability of detachment becomes small, as evident from Equation (26). In this condition the loss of the cargo becomes unlikely. Thus, the disengaged motor, on average, has enough time to attach back to the microtubule, catch up with the leading motor and move the cargo a little further leading to a net increase in the runlength. For $\overline{m}>2$, equivalent arguments can be provided and the peaks on the other curves can be similarly explained. This mechanism shows how, in the absence of the Brownian motion of the cargo, the expected runlength tends to increase while the load approaches values that are multiple of the stalling force. When the Brownian motion of the cargo is taken into account (see Figure 5 (b)) the peaks are smoothed down, but do not disappear, thus they represent a robust characteristic of the model. This kind of nonmonotonic behavior for the computed average runlength curve can be a drawback of the detachment rate model, when F is close to F_{ s }. In such a case, our approach can be seen to identify specific inaccuracies of the model. However, it is also possible that the finer analysis indicates a behavior that is exhibited by an ensemble of motors carrying a cargo and the deviation from a monotonic behavior are not artifacts of the model. In such a case, the finer analysis identifies a mechanism of “coordination” among the motors that optimizes the average runlength in situations close to a stalling scenario. Experiments can be designed and conducted in order to determine whether this mechanism is a model artifact or if it is actually occurring in the physical system. Related but different study on the relationship of velocities and runlength is reported in [23].
Detection of rare events
Thus, these “ 11 n m steps” are not associated to any actual stepping event of a motor, but exclusively to detachment events of the rearguard motor. This situation is rare in a beadassay experiment, but it is known that two motors frequently pull the shared microtubule in two opposite directions in a gliding assay experiment (see [8]). Our analysis indicates that a cargo movement with step sizes larger than 8 n m is still viable in a bead assay, though infrequent. Our approach can identify the causes for such rarer modalities of transport. Thus, steps larger that 8 n m as those described in [8] could well be originated by a mechanism of this kind. The probability distribution of the step size for an ensemble of 3 motors is reported in Figure 9 where the steps longer than 8 n m (at about 11 n m, 15 n m and 20 n m) have similar interpretation.
Conclusions
In conclusion, a framework and model for the study of the coordinated behavior of molecular motors has been introduced. The main novelty of the approach lies in the adoption of methods of analysis that obviate the need of Monte Carlo simulations.
The methodology is applied to the analysis of ensembles of KinesinI motors. Results that had been previously found using Monte Carlo methods are accurately reproduced, validating the methodology. More importantly, under this novel probabilistic approach new insights about the mechanisms of action of these proteins are found, suggesting hypothesis about their behavior and driving the design and realization of new experiments. For example, a possible mechanism under which motor proteins could coordinate together in order to increase their overall processivity is identified. Furthermore, the probabilistic framework allows the determination of steady state conditions for groups of molecular motors. The model predicts that, regardless of their initial configuration, the molecular motors will reach a situation where their relative distances on the microtubule will follow the same probability distribution. This provides an explanation for the robustness of the system with respect to the fluctuations of the surrounding environment.
The advantages provided in accuracy and efficiency make it possible to detect rare events in the motor protein dynamics, that could otherwise pass undetected using standard simulation methods. In this respect, the model has allowed to provide a possible explanation for infrequent steps of length longer that 8 n m that had been observed in bead assay experiments [8].
Declarations
Acknowledgements
The authors acknowledge NSF for its support. This research was partially supported by National Science Foundation under the grant ECCS1202411.
Authors’ Affiliations
References
 Svoboda K, Schmidt C, Schnapp B, Block S: Direct observation of kinesin stepping by optical trapping interferometry. Nature. 1993, 365 (6448): 721727. 10.1038/365721a0.View ArticleGoogle Scholar
 Milescu L, Yildiz A, Selvin P, Sachs F: Maximum likelihood estimation of molecular motor kinetics from staircase dwelltime sequences. Biophys J. 2006, 91 (4): 11561168. 10.1529/biophysj.105.079541.View ArticleGoogle Scholar
 Carter B, Vershinin M, Gross S: A Comparison of stepdetection methods: how well can you do?. Biophys J. 2008, 94: 306319. 10.1529/biophysj.107.110601.View ArticleGoogle Scholar
 Liepelt S, Lipowsky R: Operation modes of the molecular motor kinesin. Phys Rev E. 2009, 79: 011917View ArticleGoogle Scholar
 Hill D, Plaza M, Bonin K, Holzwarth G: Fast vesicle transport in PC12 neurites: velocities and forces. Eur Biophys J. 2004, 33 (7): 623632. 10.1007/s0024900404036.View ArticleGoogle Scholar
 Kural C, Kim H, Syed S, Goshima G, Gelfand VI, Selvin PR: Kinesin and dynein move a peroxisome in vivo: a tugofwar or coordinated movement?. Science. 2005, 308 (5727): 14691472. 10.1126/science.1108408.View ArticleGoogle Scholar
 Kunwar A, Vershinin M, Xu J, Gross SP: Stepping, strain gating, and an unexpected forcevelocity curve for multiplemotorbased transport. Curr Biol. 2008, 18: 117310.1016/j.cub.2008.07.027.View ArticleGoogle Scholar
 Leduc C, Ruhnow F, Howard J, Diez S: Detection of fractional steps in cargo movement by the collective operation of kinesin1 motors. Proc Nat Acad Sci. 2007, 104 (26): 1084710.1073/pnas.0701864104.View ArticleGoogle Scholar
 Soppina V, Rai AK, Ramaiya AJ, Barak P, Mallik R: Tugofwar between dissimilar teams of microtubule motors regulates transport and fission of endosomes. PNAS. 2009, 106 (46): 1938119386. 10.1073/pnas.0906524106.View ArticleGoogle Scholar
 Aggarwal T, Salapaka M: Realtime nonlinear correction of backfocalplane detection in optical tweezers. Rev Sci Instrum. 2010, 81 (12): 123105123105. 10.1063/1.3520463.View ArticleGoogle Scholar
 Meyhöfer E, Howard J: The force generated by a single kinesin molecule against an elastic load. Proc Nat Acad Sci U S A. 1995, 92 (2): 57410.1073/pnas.92.2.574.View ArticleGoogle Scholar
 Fisher M, Kolomeisky A: Simple mechanochemistry describes the dynamics of kinesin molecules. Proc Nat Acad Sci U S A. 2001, 98 (14): 774810.1073/pnas.141080498.View ArticleGoogle Scholar
 DeVille REL, VandenEijnden E: Regularity and synchrony in motor proteins. Bull Math Biol. 2008, 70 (2): 484516. 10.1007/s1153800792661.MATHView ArticleGoogle Scholar
 Klumpp S, Lipowsky R: Cooperative cargo transport by several molecular motors. Proc Nat Acad Sci U S A. 2005, 102 (48): 1728410.1073/pnas.0507363102.View ArticleGoogle Scholar
 Posta F, D’Orsogna MR, Chou T: Enhancement of cargo processivity by cooperating molecular motors. PCCP. 2009, 11: 48514860. 10.1039/b900760c.View ArticleGoogle Scholar
 Kunwar A, Mogilner A: Robust transport by multiple motors with nonlinear force–velocity relations and stochastic load sharing. Phys Biol. 2010, 7: 01601210.1088/14783975/7/1/016012.View ArticleGoogle Scholar
 Berger F, Keller C, Klumpp S, Lipowsky R: Distinct transport regimes for two elastically coupled molecular motors. Phys Rev Lett. 2012, 108 (20): 208101View ArticleGoogle Scholar
 Driver J, Rogers A, Jamison D, Das R, Kolomeisky A, Diehl M: Coupling between motor proteins determines dynamic behaviors of motor protein assemblies. Phys Chem Chem Phys. 2010, 12 (35): 1039810405. 10.1039/c0cp00117a.View ArticleGoogle Scholar
 Doob JL: Markoff chains  Denumerable case. Trans Am Math Soc. 1945, 58: 455473.MATHGoogle Scholar
 Gillespie D: Exact stochastic simulation of coupled chemical reactions. J Phys Chem. 1977, 81 (25): 23402361. 10.1021/j100540a008.View ArticleGoogle Scholar
 Leduc C, Campàs O, Zeldovich K, Roux A, Jolimaitre P, BourelBonnet L, Goud B, Joanny J, Bassereau P, Prost J: Cooperative extraction of membrane nanotubes by molecular motors. Proc Nat Acad Sci U S A. 2004, 101 (49): 1709617101. 10.1073/pnas.0406598101.View ArticleGoogle Scholar
 Beeg J, Klumpp S, Dimova R, Gracia R, Unger E, Lipowsky R: Transport of beads by several kinesin motors. Biophys J. 2008, 94 (2): 532541. 10.1529/biophysj.106.097881.View ArticleGoogle Scholar
 Xu J, Shu Z, King SJ, Gross SP: Tuning multiple motor travel via single motor velocity. Traffic. 2012, 13 (9): 11981205. 10.1111/j.16000854.2012.01385.x.View ArticleGoogle Scholar
Copyright
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.