A Stochastic Markov Model for Coordinated Molecular Motors

Donatello Materassi 1, Subhrajit Roychowdhury 1, and Murti V. Salapaka1,
1 Department of Electrical and Computer Engineering,
University of Minnesota,
200 Union St SE,
55455, Minneapolis (MN).
{mater013|subhra|murtis}@umn.edu

Abstract

Many cell functions are accomplished thanks to intracellular transport mechanisms of macromolecules along filaments. Molecular motors such as dynein or kinesin are proteins playing a primary role in these processes. The behavior of such proteins is quite well understood when there is only one of them moving a cargo particle. Indeed, numerous in vitro experiments have been performed to derive accurate models for a single molecular motor. However, in vivo macromolecules are often carried by multiple motors. The main focus of this paper is to provide an analysis of the behavior of more molecular motors interacting together in order to improve the understanding of their actual physiological behavior. Previous studies provide analyses based on results obtained from Monte Carlo simulations. Different from these studies, we derive an equipollent probabilistic model to describe the dynamics of multiple proteins coupled together and provide an exact theoretical analysis. We are capable of obtaining the probability density function of the motor protein configurations, thus enabling a deeper understanding of their behavior.

I Introduction

Many biological cell functions are accomplished thanks to intracellular transport processes. A key role in these processes is played by molecular motors that are proteins converting chemical energy into mechanical work. Among the most important molecular motors, there are kinesin and dynein. These are classes of molecular motors transporting many kind of different molecules by moving along filaments called microtubuli. They accomplish this by converting the chemical energy stored in the ATP molecules produce mechanical work and thus represent powerful nanomachines. Both kinesin and dynein have a similar structure consisting of two catalytically active motor domains, called “heads” that can bind both to a microtubule and ATP. The two heads are hinged together by a polypeptidic linkage that terminates in a distal cargo binding domain. These motor heads bind and unbind to the microtubule successively by utilizing the energy available from ATP hydrolysis and enable the motor to walk [1].

The study of molecular motors has developed over the past decade and these investigations have unraveled many previously unknown finer details at the microscopic level which have propelled various theoretical discussions of the mechanisms underlying the dynamics of molecular motors. Given the small sizes of molecular motors (each head of kinesin has a linear size of the order of 10 nm), they are supposed to be perturbed by the state of the local molecular environment and thermal fluctuations around them. But the directed walk of these cytoskeletal motors clearly reveals that they are robust enough to overcome the fluctuations surrounding them. The behavior of motor proteins is quite well understood when there is only one of them moving a cargo particle. Numerous in vitro experiments have been performed to derive accurate models to describe the motion of a single molecular motor. Indeed, it is possible to monitor the motion of a single molecular motor under highly tunable experimental conditions and obtain measurements with accurate spatial and time resolution [2, 3, 4]. As a result, there are many theoretical descriptions of motor-protein mechanisms that recognize their multiple conformational transitions taking into account the complex mechanochemical processes involved [5].
However, the chemical and mechanical properties of molecular motors in vivo are not yet completely characterized. Furthermore, there is evidence that, in vivo, macromolecules are carried by multiple motors [6]. Indeed, this can be a cause of the robustness of their behavior and of also a cause of increased efficiency. The possible dynamics of an ensemble of motor proteins carrying the same cargo particle is still largely unknown. For example, it is not clear whether synchronization phenomena occur, the motors move independently or a “tug-of-war” is in action [7, 8]. Understanding these mechanisms is one of the more challenging problems that require concerted efforts by chemists, physicists, biologists and engineers.

II Markov model for molecular motors

In this section we derive a Markov model to describe the motion of an ensemble of molecular motors carrying a cargo. While the behavior of single molecular motors such as kinesin or dynein is quite well understood, we intend to focus on the possible properties of synchronization/coordination that can be exhibited by more motors connected to the same cargo.
Molecular motors are proteins moving macro-molecules and organelles along filaments (i.e. microtubuli). Their motion occurs by discrete steps. Their heads move forward by hydrolyzing ATP and binding to proper locations that are equally spaced on the microtubule. Every motor can attach to the cargo molecule via a flexible linkage (see Figure 1).

Refer to caption — Figure 1: A molecular motor making a step

In [6] a stochastic algorithm simulating the dynamics of multiple interacting motors was employed to obtain useful insights about their behavior via Monte Carlo simulations. We “borrow” the rules of the algorithm given by [6] and introduce a rigorous mathematical formalism in order to derive a Markov model and investigate how multiple motors function together. However, in contrast to [6], our approach will have the advantage of relying on a probabilistic analysis that is exact and is not based on Monte Carlo simulations.
Every motor is modeled as a particle moving on a line by discrete steps. Its position can assume values that are equally spaced by the unit step size $d_{s}$ . Multiple motors are allowed to bind to the same location. Every motor linkage behaves as a spring with elastic constant $K$ and rest length $l_{0}$ with no compressional rigidity: when compressed it buckles down with no resistance. It is assumed that the motors are bound to the same cargo particle. The cargo is subject to the forces applied by the motors and to a constant load $F_{load}$ that is assumed positive when it opposes the motor movement.
In [6], two different models are provided for the cargo motion. In a simplified model the cargo moves deterministically reaching its equilibrium position immediately. In a more realistic model the cargo is subject to a Brownian force that models the thermal fluctuations of the environment.
In this paper we focus on the first scenario only, even though our results can be extended to treat the noisy case, as well. The motors move by single steps of length $d_{s}$ in a stochastic manner. Backward steps are not allowed. Given two non-overlapping time intervals, the probabilities that a motor takes a step in each of them are independent. This leads to a Poisson distribution of the time after which every next step occurs. For each motor, the stepping rate of the Poisson distribution depends on the reaction force $F$ it is subject to from the cargo. $F$ is assumed positive when it opposes the motor movement. A motor subject to a force $F$ will have a stepping rate given by $\lambda^{(step)}(F)$ . Estimates for $\lambda^{(step)}(F)$ can be obtained experimentally by single motor experiments and from the knowledge of the chemical process behind the ATP hydrolyzation.
Given the forces acting on each motor, they step independently from each other. The motors have a “stalling force” $F_{s}$ such that, for $F\geq F_{s}$ , we have $\lambda_{step}(F)=0$ . The stalling force can be determined via experiments, as well. The force applied to a motor also influence the probability rate $\lambda^{det}(F)$ at which it can detach form the microtubule. Every unengaged motor can also reattach to the microtubule at a rate $\lambda^{att}$ , but only at locations that are distant not more than $l_{0}$ from the cargo position. In other words, a motor can not attach to the microtubule at a location where its linkage is stretched. Finally, the model assumes that at any time there can not be more than $\overline{m}$ motors engaged on the microtubule.

II-A Mathematical formalism and Notation

We model the microtubule locations as a bi-infinite sequence $\{a_{k}\}_{k\in\mathbb{Z}}$ , $a_{k}\in\mathbb{R}$ , where each $a_{k}$ represents the position of a binding site at which a motor can be located. Since the locations are equally spaced by the distance $d_{s}$ , we provide the following definition

Definition 1

We define a location sequence as $A:=\{a_{k}\}_{k\in\mathbb{Z}}$ , where

\displaystyle a_{k}=x^{(0)}+kd_{s}

(1)

and $x^{(0)}\in\mathbb{R}$ .

Note that $x^{(0)}\in\mathbb{R}$ is determined by the arbitrary choice of the reference point on the microtubule.
The arrangement of the motors on the microtubule is described by a bi-infinite sequence of natural numbers $\{z_{k}\}_{k\in\mathbb{Z}}$ , $z_{k}\in\mathbb{N}$ , where $z_{k}$ represents the number of motors bound at the position $a_{k}$ . Equivalently, the arrangement of motors on the microtubule can be described by a finite ordered sequence $\{x_{i}\}_{i=1}^{m}$ representing the position of each of the $m$ motors engaged on the microtubule.

Definition 2

Assume a location sequence $A=\{a_{k}\}_{k\in\mathbb{Z}}$ . We refer to a configuration as any bi-infinite sequence $Z:=\{z_{k}\}_{k\in\mathbb{Z}}$ of natural numbers such that

\displaystyle m:=\sum_{k\in\mathbb{Z}}z_{k}<\infty.

(2)

The natural number $m$ is the number of engaged agents of the configuration $Z$ .
We define the position sequence of the configuration $Z$ with respect to the location sequence $A$ as the unique non-decreasing finite sequence $X:=\{x_{i}\}_{i=1}^{m}$ such that

\displaystyle z_{k}=|\{i\in\mathbb{N}\text{ such that }x_{i}\in X\text{ and }x_{i}=a_{k}\}|,

(3)

where $|\cdot|$ denotes the cardinality of a set. We write $X=\phi_{A}(Z)$ .

Each representation (configuration $Z$ or position sequence $X$ ) has notation advantages and drawbacks. For example, given a configuration, transitions like steps, attachments and detachments will be described by adding a proper “transition sequence”. Conversely, to determine the forces the cargo particle is subject to the representation it terms of a position sequence will be easier to handle.

Definition 3

We define $\mathcal{Z}_{m}$ as the set of all the possible configurations with $m$ engaged agents and $\mathcal{X}_{m}$ as the set of all possible position sequences with $m$ engaged agents. We also define

\displaystyle\mathcal{Z}:=\bigcup_{m=1}^{+\infty}\mathcal{Z}_{m}\qquad

\displaystyle\mathcal{X}:=\bigcup_{m=1}^{+\infty}\mathcal{X}_{m}.

(4)

Lemma 4

The set $\mathcal{Z}$ is countable.

Proof:

The proof is left to the reader. ∎

Lemma 5

For any $Z\in\mathcal{Z}$ , let $X=\phi_{A}(Z)\in\mathcal{X}$ be the position sequence of the configuration $Z$ with respect to the location sequence $A$ . The function $\phi_{A}:\mathcal{Z}\rightarrow\mathcal{X}$ is a one-to-one correspondence. Moreover, $\phi_{A}(\mathcal{Z}_{m})=\mathcal{X}_{m}$ for any $m\in\mathbb{N}$ .

Proof:

The proof is left to the reader. ∎

Definition 6

For any $\alpha\in\mathbb{Z}$ and $Z=\{z_{k}\}_{k\in\mathbb{Z}}\in\mathcal{Z}$ , we define the shift operator $\rho^{\alpha}:\mathcal{Z}\rightarrow\mathcal{Z}$

\displaystyle\rho^{\alpha}(Z):=\{z_{k+\alpha}\}_{k\in\mathbb{Z}}.

(5)

We assume that the motors move towards larger values of their positions. With this assumption, we introduce the following definition.

Definition 7

Consider $Z=\{z_{k}\}_{k\in\mathbb{Z}}\in\mathcal{Z}$ and let $X=\phi_{A}(Z)=(x_{1},...,x_{m})$ be its position sequence with respect to the location $A$ . We say that $x_{m}$ is the vanguard agent position while $x_{1}$ is the rearguard agent position of the configuration $Z$ (and of the position sequence $X$ ). We define spread of $Z$ (or spread of $X$ abusing the nomenclature) as

\displaystyle s(Z):=x_{m}-x_{1}

(6)

and the dimension of $Z$ as

\displaystyle h(Z):=1+\sup\{k\in\mathbb{Z}|z_{k}\neq 0\}-\inf\{k\in\mathbb{Z}|z_{k}\neq 0\}.

(7)

If $m=0$ we define $s(Z):=0$ and $h(Z):=1$ .

In the model, only certain transitions are allowed, representing a forward step, a motor detachment from the microtubule and a motor attachment to the microtubule.

Definition 8

Define the bi-infinite sequences

	$\displaystyle R^{(s)}_{k}:=\delta_{k+1}-\delta_{k}$		(8)
	$\displaystyle R^{(d)}_{k}:=-\delta_{k}$		(9)
	$\displaystyle R^{(a)}_{k}:=\delta_{k}$		(10)

where $\delta_{k}$ represents the bi-infinite Kronecker delta sequence that is equal to $1$ at $k$ and zero anywhere else.

Definition 9

Consider a time-continuous Markov Process defined on the state space $\mathcal{Z}$ , with time variable $t\in\mathbb{R}$ and with transition rate from the configuration $Z_{1}$ to $Z_{2}$ given by $\lambda(Z_{2},Z_{1})$ , meaning that, given that the system is in the configuration $Z_{1}$ at time, the probability of being in $Z_{2}$ at time $t+\Delta t$ is equal to $\lambda(Z_{2},Z_{1})\Delta t$ (for small $\Delta t$ ).
Assume that, for any $Z_{2},Z_{1}$ such that $\lambda(Z_{2},Z_{1})>0$ , there exists $k\in\mathbb{Z}$ such that one of the following conditions is met

•

$Z_{2}=Z_{1}+R^{(s)}_{k}$
•

$Z_{2}=Z_{1}+R^{(d)}_{k}$
•

$Z_{2}=Z_{1}+R^{(a)}_{k}$ .

We call the transitions sequences $R^{(s)}_{k}$ , $R^{(d)}_{k}$ and $R^{(a)}_{k}$ respectively forward step from location $k$ , detachment from location $k$ and attachment to location $k$ . We call such a process an Ensemble of Molecular Motors (EMM).

The cargo transported by the group of molecular motors is usually subject to a constant force $F_{load}$ . In this simplified article it is assumed that, after any transition, the cargo reaches its equilibrium position instantaneously. We provide the following definition.

Definition 10

Assume that a motor in position $x$ is subject to a force equal to $\chi(x-x^{(c)})$ given that the cargo is in position $x^{(c)}\in\mathbb{R}$ .
Given a position sequence $X=\{x_{1},...,x_{m}\}$ (or equivalently its configuration $Z=\phi_{A}^{-1}(X)$ ), an equilibrium position of the cargo $x^{(c)}$ is any real number satisfying the equilibrium condition

\displaystyle F_{load}=\sum_{i=1}^{m}\chi(x_{i}-x^{(c)}).

(11)

Definition 11

We refer to any realization $Z(t)$ of a EMM as a configuration path. In an obvious way, any configuration path $Z(t)$ can be associate to a position path $X(t)=\phi_{A}(Z(t))=(x_{1}(t),...,x_{m(t)}(t))$ . For any $t$ , we also have the number of engaged agents $m(t)$ , the vanguard agent with position $x_{m(t)}(t)$ , the rearguard agent $x_{1}(t)$ and, under the assumption that it is unique, the cargo equilibrium position $x^{(c)}(t)$ .

In [6] a set of specific rules are defined to create a stochastic simulation algorithm for the motion of molecular motors. The rules are the following. The cargo has $\overline{m}$ attached motors that can attach and detach from the microtubule. Their linkages behave as elastic springs of elastic constant $K$ and rest length $l_{0}$ that buckle down with no resistance when compressed. It is postulated that motors engaged on the microtubule have always a chance of detaching depending on the force they are subject to. On the other hand, any disengaged motor can attach again, but only at locations where its linkage is not subject to a stretch (that is at distance not larger than $l_{0}$ from the cargo). The probability of taking a step depends on the force acting on the specific motor and there exists a “stalling force” $F_{s}>0$ making the step not possible. Finally, the behavior of the system is invariant to translations along the microtubule. These rules are formalized in the following way

Definition 12

Given a EMM, we say that it satisfies the Kunwar-Veshinin-Xu-Gross (KVXG) rules with parameters $(\overline{m},K,l_{0},F_{s})$ , with $\overline{m}\in\mathbb{N}$ and $K,l_{0},F_{s}$ real positive, if it meets the following properties. Assume a location sequence $A=\{a_{k}\}_{k\in Z}$ ; for any $Z_{1}\neq Z_{2}\in\mathcal{Z}$ , with position sequences $X_{1}$ , $X_{2}$ and engaged agents $m_{1}$ , $m_{2}$ respectively, we have

•

irreversible cargo loss with no engaged motors

$\displaystyle m_{1}=0\text{ implies }\lambda(Z_{2},Z_{1})=0$ (12)
•

at most $\overline{m}$ engaged motors

$\displaystyle m_{2}>\overline{m}\text{ implies }\lambda(Z_{2},Z_{1})=0$ (13)
•

spatial invariance (for any $\alpha\in\mathbb{Z}$ )

$\displaystyle\lambda(Z_{2},Z_{1})=\lambda(\rho^{\alpha}Z_{2},\rho^{\alpha}Z_{1})$ (14)

•

elastic stretch with constant $K$ and dead zone $l_{0}$

\displaystyle\chi(l):=\left\{\begin{array}[]{ll}K(l-l_{0})&\text{if }l>l_{0}\\ 0&\text{if }|l|\leq l_{0}\\ K(l+l_{0})&\text{if }l<-l_{0}.\end{array}\right.

(18)

•

attachment with no stretch

	$\displaystyle\lambda(Z_{2},Z_{1})>0\text{ and }Z_{2}-Z_{1}=R^{(a)}_{k}$
	$\displaystyle\text{ implies }\|a_{k}-x^{(c)}\|\leq l_{0}$		(19)

•

stalling force $F_{s}$

$\displaystyle\chi(a_{k}-x^{(c)})\geq F_{s}\text{ implies }\lambda(Z_{1}+R^{(s)}_{k},Z_{1})=0.$ (20)

Given an EMM with a transition rate function $\lambda:\mathcal{Z}\rightarrow\mathbb{R}^{+}$ meeting the KVXG rules, define $P_{Z}(Z,t|\overline{Z},t_{0})$ as the probability of being in the configuration $Z=\{z_{k}\}_{k\in\mathbb{Z}}$ at time $t\geq t_{0}$ conditioned to $Z(t_{0})=\overline{Z}$ .

Then, given an initial time $t_{0}$ and an initial state $\overline{Z}$ , for $t\geq t_{0}$ , $P_{Z}(Z,t|\overline{Z},t_{0})$ satisfies the Master Equation

		$\displaystyle\frac{\partial}{\partial t}P_{Z}(Z,t\|\overline{Z},t_{0})=-P_{Z}(Z,t\|\overline{Z},t_{0})\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(Z^{\prime},Z)$		(21)
		$\displaystyle\qquad+\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(Z,Z^{\prime})P_{Z}(Z^{\prime},t\|\overline{Z},t_{0}),$		(22)

that represents the conservation law of the probability measure.

III Finite dimensional projection and Master Equation analysis

The study of the dynamics of the motors (represented either in the configuration space or in the position sequence space) is not easy to handle since the entities involved are not stationary and the state-space is infinite-dimensional. The following theoretical result is helpful for transforming the problem into a finite dimensional one.

Theorem 13

Consider an EMM satisfying the KVXG rules with parameters $(\overline{m},K,l_{0},F_{s})$ . For any configuration path $Z(t)$ let $s(Z(t))$ be its associated spread for any $t$ . Define

\displaystyle s^{(max)}:=\max\left\{\frac{\overline{m}F_{s}-F_{load}}{K}+d_{s},\frac{F_{load}}{K}\right\}+2l_{0}.

(23)

For any $S\geq s^{(max)}$ and for any $t\geq t_{0}$ , we have that if $s(\overline{Z})\leq S$

\displaystyle Pr\{s(Z(t))\leq S~|~Z(t_{0})=\overline{Z}\}=1.

(24)

where $Pr\{E|C\}$ is the conditional probability of an event $E$ given the condition $C$ .

Proof:

See the appendix ∎

Theorem 13 enables us to study many properties of an EMM using a formulation that has finite dimension and deals with stationary quantities. To this aim, we introduce the following definitions.

Definition 14

Given an EMM satisfying the KVXG rules with parameters $(\overline{m},K,l_{0},F_{s})$ subject to load $F_{load}$ and a location sequence $A$ with step size $d_{s}$ , we define

\displaystyle n:=\left\lceil\frac{s^{(max)}}{d_{s}}\right\rceil.

(25)

as the regular dimension of the ensemble.
Given a configuration $Z=\{z_{k}\}_{k\in\mathbb{Z}}$ , we also define its ensemble representation as the $n$ -vector

\displaystyle Q=(z_{k_{1}},...z_{k_{1}+n-1})^{T}

(26)

where $k_{1}:=\inf\{k\in\mathbb{Z}|z_{k}\neq 0\}$ and we write $Q:=\Pi^{(e)}(Z)$ . If $k_{1}=-\infty$ , $Q$ is $(0,0,...,0)^{T}$ .

The ensemble representation of $Z$ is the projection of the configuration on a suitable $n$ -dimensional space. Theorem 13 guarantees that the ensemble representation $Q(t)$ of the configuration $Z(t)$ contains all the non-zero entries of $Z(t)$ with probability $1$ for all $t>t_{0}$ if the dimension of $Z(t_{0})$ is not greater than $n$ . In this scenario, while $Q(t)$ does not contain the information about the absolute positions of the motors, it still contains all the information about all their relative positions, so it can be used as a tool to study the dynamics of their mutual interactions.
Assume that the starting configuration $Z(t_{0})=\overline{Z}$ has dimension not exceeding the regular one $n$ . For $t>t_{0}$ the spread is not exceeding $n$ and the transition probabilities of $Z(t)$ are known. Despite the projection, it is possible to derive the transition probabilities of $Q(t)$ . This allows to define the dynamics of a Markov system with a finite dimension. The following properties hold.

Lemma 15

Consider an EMM with regular dimension $n$ . Let $s(\overline{Z})\leq n$ and $P_{Q}(Q,t)$ be the probability that $\Pi^{(e)}(Z)$ is equal to $Q$ at time $t$ . Then, for $Q\neq 0$ and for any $t\geq t_{0}$

\displaystyle P_{Q}(Q,t|\overline{Z},t_{0})=\sum_{\alpha\in\mathbb{Z}}P_{Z}\{\rho^{\alpha}Z,t|\overline{Z},t_{0}\}

(27)

Proof:

The proof is left to the reader. ∎

Theorem 16

Consider an EMM with regular dimension $n$ and let $s(t)$ be the spread at time $t$ . Define

\displaystyle\lambda_{Q}(Q^{\prime},Q):=\sum_{\stackrel{{\scriptstyle\Pi^{(e)}(Z^{\prime})=Q^{\prime}}}{{s(Z)\leq n}}}\lambda(Z^{\prime},Z).

(28)

For any $t\geq t_{0}$ and for any $\overline{Z}$ such that $s(\overline{Z})\leq n$ , it holds that

	$\displaystyle\frac{\partial}{\partial t}P_{Q}(Q,t\|\overline{Z},t_{0})=-P_{Q}(Q,t\|\overline{Z},t_{0})\sum_{Q^{\prime}\in\mathcal{Q}}\lambda(Q^{\prime},Q)$		(29)
	$\displaystyle\qquad+\sum_{Q^{\prime}\in\mathcal{Q}}\lambda(Q,Q^{\prime})P_{Q}(Q^{\prime},t\|\overline{Z},t_{0})$		(30)

Proof:

In order to simplify the notation we omit the initial condition: $P_{Q}(Q,t)=P_{Q}(Q,t|\overline{Z},t_{0})$ and $P_{Z}(Z,t)=P_{Z}(Z,t|\overline{Z},t_{0})$ . We have that $Q=0$ if and only if $Z=0$ . In such a case, $\lambda(Z^{\prime},Z)=0$ for any $Z^{\prime}$ and $\lambda(Z,Z^{\prime})=0$ for any $Z^{\prime}$ such that $\Pi^{(e)}(Z^{\prime})\neq Q_{1}:=(1,0,...,0)$ . Observe that

•

$\lambda_{Q}(Q^{\prime},0)=0$ for any $Q^{\prime}$
•

$\lambda_{Q}(0,Q^{\prime})=0$ for any $Q^{\prime}\neq Q_{1}$
•

$P_{Z}(Z^{\prime},t)=0$ if $s(Z^{\prime})>n$ .

Then, for $Q=0$ we obtain

	$\displaystyle\frac{\partial}{\partial t}P_{Q}(0,t)=\frac{\partial}{\partial t}P_{Z}(0,t)=$
	$\displaystyle=\sum_{\stackrel{{\scriptstyle\Pi^{(e)}(Z^{\prime})=Q_{1}}}{{s(Z)\leq n}}}\lambda(0,Z^{\prime})P_{Z}(Z^{\prime},t)=\lambda_{Q}(0,Q_{1})P_{Q}(0,Q_{1}).$

For $Q\neq 0$ , we have

	$\displaystyle\frac{\partial}{\partial t}P_{Q}(Q,t)=\frac{\partial}{\partial t}P_{Q}(\Pi^{(e)}(Z),t)=\sum_{\alpha\in\mathbb{Z}}\frac{\partial}{\partial t}P_{Z}(\rho^{\alpha}Z,t)=$
	$\displaystyle\qquad-\sum_{\alpha\in\mathbb{Z}}\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(Z^{\prime},\rho^{\alpha}Z)P_{Z}(\rho^{\alpha}Z,t)$
	$\displaystyle\qquad+\sum_{\alpha\in\mathbb{Z}}\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(\rho^{\alpha}Z,Z^{\prime})P_{Z}(Z^{\prime},t)=$
	$\displaystyle=-\sum_{\alpha\in\mathbb{Z}}\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(\rho^{\alpha}Z^{\prime},\rho^{\alpha}Z)P_{Z}(\rho^{\alpha}Z,t)$
	$\displaystyle\qquad+\sum_{\alpha\in\mathbb{Z}}\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(\rho^{\alpha}Z,\rho^{\alpha}Z^{\prime})P_{Z}(\rho^{\alpha}Z^{\prime},t)$
	$\displaystyle=-\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(Z^{\prime},Z)\sum_{\alpha\in\mathbb{Z}}P_{Z}(\rho^{\alpha}Z,t)$
	$\displaystyle\qquad+\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(Z,Z^{\prime})\sum_{\alpha\in\mathbb{Z}}P_{Z}(\rho^{\alpha}Z^{\prime},t)=$
	$\displaystyle=-P_{Q}(Q,t)\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(Z^{\prime},Z)$
	$\displaystyle\qquad+\sum_{\stackrel{{\scriptstyle Z^{\prime}\in\mathcal{Z}}}{{s(Z^{\prime})\leq n}}}\lambda(Z,Z^{\prime})P_{Q}(\Pi^{(e)}(Z^{\prime}))=$
	$\displaystyle=-P_{Q}(Q,t)\sum_{Q^{\prime}\in\mathcal{Q}}\lambda_{Q}(Q^{\prime},Q)$
	$\displaystyle\qquad+\sum_{Q^{\prime}\in\mathcal{Q}}\lambda_{Q}(Q,Q^{\prime})P_{Q}(Q^{\prime})$

∎

Lemma 17

Given an EMM satisfying the KVXG rules with parameters $(\overline{m},K,l_{0},F_{s})$ with regular dimension $n$ , the number of possible ensemble representations is given by

\displaystyle N=1+\sum_{m=1}^{\overline{m}}\frac{(n+m-2)!}{(n-1)!(m-1)!}.

(31)

Proof:

The proof is left to the reader. ∎

Enumerate all the possible ensemble representations $Q_{1},...,Q_{N}$ with $N$ given by (31) and define $P_{i}(t)$ as the probability of having the ensemble representation $Q_{i}$ at time $t$ given the initial condition $\overline{Z}$ such that $s(\overline{Z})\leq n$ . By projecting Equation (21) onto the set $\mathcal{Q}:=\{Q_{1},...,Q_{N}\}$ we obtain

	$\displaystyle\frac{d}{dt}P_{Q}(Q,t)=-\sum_{Q^{\prime}\in\mathcal{Q}^{(+)}(Q)}\lambda_{Q}(Q^{\prime},Q)P_{Q}(Q,t)$		(32)
	$\displaystyle+\sum_{Q^{\prime}\in\mathcal{Q}^{(-)}(Q)}\lambda_{Q}(Q,Q^{\prime})P_{Q}(Q^{\prime},t),$		(33)

The dynamics of the probabilities on the space of the ensemble representations is given by an equation of the form

\displaystyle\frac{d}{dt}P=AP

(34)

where $P:=(P_{1},...,P_{N})^{T}$ , the off-diagonal element of $A$ , $a_{ji}$ is the rate of the transition from $Q_{i}$ to $Q_{j}$ and the sum of all the rows of $A$ is zero. In this way the dynamics of the infinite dimensional Markov system has been projected onto the dynamics of a finite dimensional one. In the section we show how this projection provides no loss of information when used to obtain certain quantities of interest such as the expected velocity of the cargo and the expected runlength.

IV Determining the transition probabilities from single motor observations

In this section we determine the transition rates $\lambda(Z^{\prime},Z)$ for the infinite dimensional stochastic model. The transition rates $\lambda_{Q}(Q^{\prime},Q)$ necessary to define the finite dimensional one can be computed from the transition rates $\lambda(Z^{\prime},Z)$ . Following [6], it is possible to derive the probability rates of simple events $\lambda(Z^{\prime},Z)$ by considering the behavior of a single motor.

Rate for stepping transitions

Every motor moves converting ATP into kinetic energy according to the chemical equation

K+ATP\stackrel{{\scriptstyle[}}{{k}}_{off}]{k_{on}}{\rightleftarrows}K~ATP\stackrel{{\scriptstyle k_{cat}}}{{\rightarrow}}K+ADP+P_{i}+{\tt energy}.

(35)

Following [9], a Michaelis-Menten dynamics leads to a hydrolisis rate equal to $k_{cat}[ATP]/([ATP]+k_{m})$ . It is also assumed that the free head of the motor successfully binds to the microtubule location with probability rate (or efficiency) $\varepsilon$ . Then, the transition rate $\lambda_{step}$ of stepping and the average velocity $V$ of a single motor are given by

	$\displaystyle\lambda_{step}=\frac{k_{cat}[ATP]}{[ATP]+k_{m}}\varepsilon$		(36)
	$\displaystyle V=P_{step}d=\frac{k_{cat}[ATP]}{[ATP]+k_{m}}d\varepsilon$		(37)

where $d$ is the fixed size of the protein step. The force $F$ that the cargo exerts on the motor is assumed positive when it opposes the motor motion. It is assumed that when it exceeds a certain force $F_{0}$ it causes the motor to stall. Following [6], it is assumed to affect the motor dynamics by changing the probability $\varepsilon$ of binding with the microtubule, following the relation

\displaystyle\varepsilon(F)=\left\{\begin{array}[]{lr}1&\text{if }F\leq 0\\ 1-\left(\frac{F}{F_{s}}\right)^{2}&\text{if }0<F<F_{0}\\ 0&\text{otherwise}\end{array}\right.

(41)

where $F_{s}$ is the stalling force of the KVXG rules. In Figure 2, it is represented how the efficiency $\varepsilon(F)$ changes in different situations: when the motor is lagging behind the cargo, when it is at a distance where the linkage is not stretched and when it is actively pulling the cargo.

The linkage is modeled as an elastic spring of elastic constant $K$ and rest length $l_{0}$ that buckles down with no resistance when compressed. The relation between the force applied to the motor and the linkage length $l$ is given by

\displaystyle F(l)=\left\{\begin{array}[]{lr}K(l+l_{0})&\text{if }l\leq-l_{0}\\ 0&\text{if }|l|<l_{0}\\ K(l-l_{0})&\text{if }l\geq l_{0}.\\ \end{array}\right.

(45)

In [6] it is assumed that the force $F$ also influences the kinetics of the ATP hydrolysis. In particular it is assumed that $k_{off}$ increases with increasing $F$

k_{off}=k_{0off}e^{Fd_{l}/K_{b}T},

(46)

where $k_{0off}$ is the backward reaction rate of the hydrolysis when $F=0$ , $K_{b}$ is the Boltzmann constant, $T$ is the temperature and $d_{l}$ is a parameter that can be experimentally determined. Thus, the transition rate for a step at location $a_{k}$ is given by

		$\displaystyle\lambda(Z+R_{k}^{(s)},Z)=z_{k}\frac{k_{cat}[ATP]}{[ATP]+\frac{k_{on}+k_{off}(F)}{k_{cat}}}\epsilon(F)$		(47)
		$\displaystyle=z_{k}\frac{k_{cat}[ATP]\left[1-\left(\frac{\chi(a_{k}-x^{(c)})}{F_{0}}\right)^{2}\right]}{[ATP]+\frac{k_{on}+k_{0off}e^{\chi(a_{k}-x^{(c)})d_{l}/K_{b}T}}{k_{cat}}}$		(48)

Rate for detachment transitions

From [Schnitzer et al. 2000], the processivity $L$ is

\displaystyle L=\frac{d[ATP]Ae^{-F\delta_{l}K_{b}T}}{[ATP]+B(1+A)e^{-F\delta_{l}K_{b}T}},

(49)

where $A$ , $B$ and $\delta_{l}$ are parameters that can be experimentally determined. Since the processivity represents how far a motor can move before detaching from the microtubule on average, we find a relation between the probability of stepping and the probability of detachment.

\displaystyle\frac{P_{step}(F)}{P_{detach}(F)}=\frac{L}{d}=\frac{[ATP]Ae^{-F\delta_{l}K_{b}T}}{[ATP]+B(1+A)e^{-F\delta_{l}K_{b}T}}.

(50)

Thus, so long as $F<F_{0}$ ,

\displaystyle P_{detach}(F)=\frac{[ATP]+B(1+A)e^{-F\delta_{l}K_{b}T}}{[ATP]Ae^{-F\delta_{l}K_{b}T}}P_{step}(F).

(51)

When $F\geq F_{s}$ , in [6] a constant detachment rate is assumed $P_{detach}(F)=P_{back}=2s^{-1}$ . Thus, we have the following transition rates

\displaystyle\lambda(Z+R_{k}^{(d)},Z)=\left\{\begin{array}[]{l}\frac{[ATP]+B(1+A)e^{-\chi(a_{k}-x^{(c)})\delta_{l}K_{b}T}}{[ATP]Ae^{-\chi(a_{k}-x^{(c)}\delta_{l}K_{b}T}}\lambda(Z+R_{k}^{(s)},Z)\\ \qquad\text{ if }\chi(a_{k}-x^{(c)}<F_{s}\\ z_{k}P_{back}\\ \qquad\text{ if }\chi(a_{k}-x^{(c)}\geq F_{s}\end{array}\right.

(56)

Probability of attachment

Experimentally, it is found in [6] that the probability of a motor to attach to the microtubule is $P_{att}\simeq 5s^{-1}$ . If the motor is linked to the cargo, it is assumed that is 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.

V Probabilistic Analysis

We have shown that the length of all the possible strings can be uniformly bounded. As a consequence, apart from possible spatial shifts along the microtubule, the number of possible arrangements of the motors is finite.

Using the finite dimensional Markov Model derived in the previous sections, it is possible to derive meaningful characteristics of the system in an exact manner.

V-A Transient analysis

In Figure 3 we present the dynamics of $P(t)$ in the case $\overline{m}=3$ for two different values of $F_{load}$ . We use $4$ time snapshots to describe how $P(t)$ changes with time.

In the first row of Figure 3 we present the probability density function of the configurations when the cargo is subject to a low load ( $F_{load}=0.0002nN$ ). The $x$ and the $y$ axis represent the distance of the two leading motors from the rearguard motor. The distance is expressed in number of locations on the microtubule. The $z$ axis represents the probability of having a motor at distance $x$ and a motor a distance $y$ from the rearguard motor. The probabilities of configurations with less than three motors engaged are not reported because that would make the visualization unnecessarily more complicated. The four snapshots are taken at time $0sec$ , $0.2sec$ , $0.4sec$ , $0.6sec$ from left to right. In the second row of Figure 3 we report the results in the case of a high load ( $F_{load}=0.008nN$ ) with an identical initial probability density function. Notice how the change in the parameter $F_{load}$ leads to two different dynamics of the probability density function: under a low load the motors tend to spread out while under a high load certain configurations tend to be more likely.

V-A1 Average runlength

The detailed knowledge of the probability vector $P(t)$ can be used to derive information about numerous aspects of the behavior of the motor ensemble and study the variation of certain quantities as a function of different parameters. An important quantity that can also be experimentally measured in experiments is the expected runlength of the motors, that is the average length travelled by the cargo before being lost. The average runlength is immediately computed from $P(t)$

\displaystyle Average~Runlength=E\left[\int_{0}^{+\infty}\frac{dx^{(c)}}{dt}dt\right]

(57)

where $E[\cdot]$ is the mean operator and $x^{(c)}$ is the cargo position. This same quantity as a function of the load has already been computed in [6] using Monte Carlo methods in two scenarios. In one scenario (Model A) the authors neglect the effect of thermal noise and in another scenario (Model B) they provide a dynamic model for the brownian motion of the cargo.
We report our results in Figure 4 for two reasons. The first reason is for comparison. By forcing the variance of the cargo position to zero we obtain a model that is equivalent to “Model A” in [6]. The only difference is that in our case we are computing the average runlength from the probability density function without using a simulation approach. Thus, forcing the variance of the cargo position to zero, the results must match. We report the results of our computation in Figure 4 a using both a coarse grid (solid lines) and a fine grid (dashed lines) for the load force. If we compare the results we computed with the results in [6] at the same points (that is using the coarse grid), the runlength-load curves overlap almost perfectly.


(a)	(b)

The second reason is to highlight a possible weakness of the model adopted for the detachment rate. In the finer scale, we can notice peaks in the runlength curves that were not evident before. These peaks appear in correspondence of loads that are multiple of the stalling force $F_{0}$ and can be explained as “artifacts” of the model. They can be explained in the following way. Let us consider the case $\overline{m}=2$ for simplicity. When only one motor is engaged ( $m(t)=1$ ) and $F_{load}$ is close to (but less than) the stalling force $F_{0}$ , the probability of detachment becomes small, and the loss of the cargo becomes unlikely. In the meantime, the disengaged motor has the opportunity of attaching to microtubule, catching up with the leading motor and moving the cargo a little further. For $\overline{m}>2$ , similar arguments can be used. This mechanisms explains how, in the mathematical model, by neglecting the brownian motion of the cargo, the expected runlength tends to increase while the load approaches multiple values of the stalling force.
When the brownian motion of the cargo is taken into account (see Figure 4 b) those peaks are smoothed down, but do not disappear completely. Thus, the pronounced peaks in Figure 4 a and Figure 4 b are due to the inaccuracy of the equation (47) for $F\simeq F_{s}$ . Still, the insights given by the model could have a physical meaning and could be used to explain possible mild breaks of monotonicity in the runlength curves for multiple motors.

V-B Steady state characteristics

V-B1 Average number of engaged motors

In Figure 5 a, the probability distribution $P(t)$ has been used to determine the probability of having a given number of motors engaged on the microtubule at time $t$ .

We notice that the probability of cargo loss increases with time (and it can be proved converges to $1$ for $t\rightarrow+\infty$ ). Thus, it is not possible to define a steady state for the system: the eventual cargo loss is the actual steady state. However, in Figure 5 b we report the conditional probability density function of the number of motors engaged on the microtubule at time $t$ , given that at least one motor is engaged. We observe that using this conditional probability density function a steady state can be defined. In the case of $\overline{m}=3$ , the steady state conditional probability is reported in the histograms of Figure 6 a and in Figure 6 b for $F{load}=0.0002nN$ and $F{load}=0.008nN$ respectively. As in Figure 3, the rearguard motor is taken as a reference point, the $x$ and $y$ axis represent the distance of the other two motors from the rearguard motor and the $z$ axis represents the probability.

It can be observed that under a low load the motors tend to be spread and distant from each other, while under a high load they tend to be clusterized. Indeed, the pronounced probability peak around the point $x=0,y=0$ indicates that the three motors are not far from each other. It is interesting to notice that there is a significantly higher probability along the diagonal line $x=y$ . This can be interpreted as follows. The three motors tend to be close together, but when they are not it is more likely to find the two leading motors together while the rearguard motor is lagging behind.

V-B2 Average velocity

The definition of the steady state conditional probability can be used to determine the average steady state velocity of the motor ensemble. The results are reported in Figure 7 a for the case where thermal fluctuations are neglected and in Figure 7 b for the more accurate model. Results in Figure 7 a are again for comparison and verification since, in the noiseless scenario, our model is equivalent to Model A in [6]. Results qualitatively and quantitatively non-dissimilar are reported in Figure 7 b for the noisy scenario.

V-C Detection of rare events

The possibility of determining exactly the probability distribution of the different motor configurations allows for the detection of rare events. It is possible, for example, to determine the probability of the different steps sizes for an ensemble of $2$ motors as represented in Figure 8.

As it can be noted, there are two prominent peaks corresponding to $8~nm$ and $4~nm$ . These peaks corresponds respectively to the case where there is one active motor before and after the step and to the case where there are two active motors before and after the step. There are also different predicted step sizes at about $2~nm$ , $6~nm$ and around $4~nm$ . They correspond to situations where there is a different number of active motors before and after the step and they can be mostly considered as artifacts of the model. However, we also find a small probability, but different from zero, of steps larger than $8~nm$ at about $11~nm$ . This does not seem possible, since the each motor can advance only $8~nm$ at the time. Indeed, this can occur when the rearguard motor is pulling back the cargo and detaches from the microtubule. While this situation is not very likely in a bead-assay experiment that our simulations are trying to portray, it is quite likely that two motors could pull in two opposite directions in a gliding assay experiment. Thus, steps larger that $8~nm$ as those described in [10] could 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 8 where the steps longer than $8~nm$ have an identical interpretation.

-D Preliminary Lemmas

Lemma 18

If $F_{load}>0$ and $m>0$ , then the equilibrium solution $x^{(c)}$ to (11) for any position sequence $X=(x_{1},...,x_{m})$ and for $\chi$ as given by (18) exists, is unique and $x_{m}-l_{0}>x^{(c)}$ .

Proof:

The proof is left to the reader. ∎

Lemma 19

Consider an EMM meeting the KVXG rules with parameters $(\overline{m},K,l_{0},F_{s})$ . Assume that $F_{load}>0$ , and that, for a configuration $Z$ , the number of engaged agents is $m>0$ such that the cargo position $x^{(c)}$ is uniquely defined. Also assume $X=\phi_{A}(Z)=(x_{1},...,x_{m})$ for a location sequence $A$ and that

s(Z):=x_{m}-x_{1}\geq\frac{mF_{s}-F_{load}}{K}+2l_{0}.

(58)

Then the vanguard motor is stalled, that is $F_{s}\leq\chi(x_{m}-x^{(c)})$ . If the vanguard motor is in the location $a_{k_{m}}$ it immediately follows that, for $k\geq k_{m}$ , $\lambda(Z+R^{(s)}(k),Z)=0$ .

Proof:

Without any loss of generality, assume $x^{(c)}=0$ , since the reference point on the microtubule is arbitrary. By contradiction, and by using Lemma (18), assume $K(x_{m}-l_{0})<F_{s}$ . If $x_{1}<-l_{0}$ , then we have

	$\displaystyle F_{load}-K(x_{1}+l_{0})=F_{load}+\sum_{i=2}^{m}\chi(x_{i})\leq$		(59)
	$\displaystyle\leq(m-1)K(x_{m}-l_{0})<-K(x_{m}-l_{0})+mF_{s}$		(60)

This implies

\displaystyle x_{m}-x_{1}<\frac{mF_{s}-F_{load}}{K}+2l_{0}

(61)

that is a contradiction. In the case $x_{1}\geq-l_{0}$ , we have

\displaystyle F_{load}=\sum_{i=1}^{m}\chi(x_{i})\leq mF_{s}<mF_{s}+K(x_{m}-l_{0}).

(62)

implying

\displaystyle x_{m}-x_{1}\leq x_{m}+l_{0}<\frac{mF_{s}-F_{load}}{K}+2l_{0},

(63)

that is again a contradiction. ∎

Lemma 20

\displaystyle q\geq\frac{F_{load}}{K}+2l_{0},

(64)

and that the next transition occurring is an attachment at position $k$ , that is $Z^{\prime}=Z+R^{(a)}(k)$ . Then, if $\lambda(Z^{\prime},Z)>0$

\displaystyle s(Z^{\prime})\leq q.

(65)

Proof:

Without any loss of generality, assume $x^{(c)}=0$ , since the reference point on the microtubule is arbitrary. Let $X^{\prime}=(x^{\prime}_{1},...,x_{m+1}^{\prime})=\phi_{A}(Z^{\prime})$ If $x_{1}<-l_{0}$ , the transition does not affect the value of the distance because the newly attached motor will be in the range $[-l_{0},l_{0}]$ in order to meet the KVXG rules. Indeed, we have

\displaystyle s(Z^{\prime})=x_{m+1}^{\prime}-x_{1}^{\prime}=x_{m}-x_{1}=s(Z)\leq q

(66)

If $x_{1}\geq-l_{0}$ , then $x_{1}^{\prime}\geq-l_{0}$ . Moreover,

\displaystyle F_{load}=\sum_{i=1}^{m}\chi(x_{i})\geq K(x_{m}-l_{0}),

(67)

implies

\displaystyle x_{m+1}^{\prime}-x_{1}^{\prime}\leq x_{m}+l_{0}\leq\frac{F_{load}}{K}+2l_{0}\leq q.

(68)

∎

-E Proof of Theorem 13

Proof:

Define, for $t\geq t_{0}$ ,

	$\displaystyle 0\leq\phi(t):=Pr\{s(t)>S\|Z(t_{0})=\overline{Z}\}=$		(69)
	$\displaystyle\qquad=\sum_{s(Z)>S}P_{Z}(Z,t\|\overline{Z},t_{0}).$		(70)

At time $t_{0}$ we have that $\phi(t_{0})=0$ . By contradiction assume that there exists a time $t_{1}>t_{0}$ such that $\phi(t_{1})>0$ . We must also have that, for some $t\in(t_{0},t_{1})$ ,

\displaystyle\frac{\partial}{\partial t}\phi(t)>0.

(71)

From the master equation

	$\displaystyle\frac{\partial}{\partial t}\phi(t)=\sum_{s(Z)>S}\left\{-\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(Z^{\prime},Z)P_{Z}(Z,t\|\overline{Z},t_{0})\right.+$
	$\displaystyle\qquad\left.+\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(Z,Z^{\prime})P_{Z}(Z^{\prime},t\|\overline{Z},t_{0})\right\}=$
	$\displaystyle=\sum_{s(Z)>S}\sum_{s(Z^{\prime})\leq S}\left\{-\lambda(Z^{\prime},Z)P_{Z}(Z,t\|\overline{Z},t_{0})\right.+$
	$\displaystyle\qquad\left.+\lambda(Z,Z^{\prime})P_{Z}(Z^{\prime},t\|\overline{Z},t_{0})\right\}\leq$
	$\displaystyle\leq\sum_{s(Z)>S}\sum_{s(Z^{\prime})\leq S}\left\{\lambda(Z,Z^{\prime})P_{Z}(Z^{\prime},t\|\overline{Z},t_{0})\right\}$

Now we prove that, if $s(Z^{\prime})\leq S$ and $s(Z)>S$ , we necessarily have $\lambda(Z,Z^{\prime})=0$ .
The only possible transitions are detachment, attachment and step.
The detachment transition can only reduce the configuration distance, so $Z-Z^{\prime}\neq R^{(d)}(k)$ for any $k$ .
If the transition is an attachment, by Lemma 20 we obtain the assertion.
If the transition is a step the configuration distance increases only if the vanguard motor steps, and the increase is by $d_{s}$ . If $s(Z^{\prime})<S-d_{s}$ , then $s(Z)<S$ and this is not a case we are considering because we have assumed $s(Z)>S$ . Conversely, if $s(Z^{\prime})\geq S-d_{s}$ , by Lemma 19, we have that the vanguard motor can not step, so such a transition is not possible. Then we obtain as a contradiction that, for $t\geq t_{0}$ ,

\displaystyle\frac{\partial}{\partial t}\phi(t)\leq 0.

(72)

∎

References

[1] N. J. Carter and R. Cross, “Mechanics of the kinesin step,” Nature, vol. 435, no. 7040, p. 308, 2005.
[2] K. Svoboda, C. Schmidt, B. Schnapp, and S. Block, “Direct observation of kinesin stepping by optical trapping interferometry,” Nature, vol. 365, no. 6448, pp. 721–727, 1993.
[3] L. Milescu, A. Yildiz, P. Selvin, and F. Sachs, “Maximum likelihood estimation of molecular motor kinetics from staircase dwell-time sequences,” Biophysical journal, vol. 91, no. 4, pp. 1156–1168, 2006.
[4] B. Carter, M. Vershinin, and S. Gross, “A Comparison of Step-Detection Methods: How Well Can You Do?” Biophysical Journal, vol. 94, no. 1, pp. 306–319, 2008.
[5] S. Liepelt and R. Lipowsky, “Operation modes of the molecular motor kinesin,” Phys. Rev. E, vol. 79, p. 011917, 2009.
[6] A. Kunwar, M. Vershinin, J. Xu, and S. P. Gross, “Stepping, strain gating, and an unexpected force-velocity curve for multiple-motor-based transport,” Current biology, vol. 18, p. 1173, 2008.
[7] C. Kural, H. Kim, S. Syed, G. Goshima, V. I. Gelfand, and P. R. Selvin, “Kinesin and dynein move a peroxisome in vivo: a tug-of-war or coordinated movement?” Science, vol. 308, no. 5727, pp. 1469–1472, 2005.
[8] V. Soppina, A. K. Rai, A. J. Ramaiya, P. Barak, and R. Mallik, “Tug-of-war between dissimilar teams of microtubule motors regulates transport and fission of endosomes,” PNAS, vol. 106, no. 46, pp. 19 381–19 386, 2009.
[9] E. Meyhöfer and J. Howard, “The force generated by a single kinesin molecule against an elastic load,” Proceedings of the National Academy of Sciences of the United States of America, vol. 92, no. 2, p. 574, 1995.
[10] C. Leduc, F. Ruhnow, J. Howard, and S. Diez, “Detection of fractional steps in cargo movement by the collective operation of kinesin-1 motors,” Proceedings of the National Academy of Sciences, vol. 104, no. 26, p. 10847, 2007.

	$\displaystyle\frac{\partial}{\partial t}\phi(t)=\sum_{s(Z)>S}\left\{-\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(Z^{\prime},Z)P_{Z}(Z,t\|\overline{Z},t_{0})\right.+$
	$\displaystyle\qquad\left.+\sum_{Z^{\prime}\in\mathcal{Z}}\lambda(Z,Z^{\prime})P_{Z}(Z^{\prime},t\|\overline{Z},t_{0})\right\}=$
	$\displaystyle=\sum_{s(Z)>S}\sum_{s(Z^{\prime})\leq S}\left\{-\lambda(Z^{\prime},Z)P_{Z}(Z,t\|\overline{Z},t_{0})\right.+$
	$\displaystyle\qquad\left.+\lambda(Z,Z^{\prime})P_{Z}(Z^{\prime},t\|\overline{Z},t_{0})\right\}\leq$
	$\displaystyle\leq\sum_{s(Z)>S}\sum_{s(Z^{\prime})\leq S}\left\{\lambda(Z,Z^{\prime})P_{Z}(Z^{\prime},t\|\overline{Z},t_{0})\right\}$