Abstract

A Deterministic Model for One-Dimensional Excluded Flow with Local Interactions

Yoram Zarai, Michael Margaliot, Anatoly B. Kolomeisky

1 School of Electrical Engineering, Tel-Aviv University, Tel-Aviv 69978, Israel.

2 School of Electrical Engineering and the Sagol School of Neuroscience, Tel-Aviv University, Tel-Aviv 69978, Israel.

3 Department of Chemistry, Rice University, Houston, TX 77005-1892, USA.

Corresponding Author. Email: michaelm@post.tau.ac.il

## Abstract

Natural phenomena frequently involve a very large number of interacting molecules moving in confined regions of space. Cellular transport by motor proteins is an example of such collective behavior. We derive a deterministic compartmental model for the unidirectional flow of particles along a one-dimensional lattice of sites with nearest-neighbor interactions between the particles. The flow between consecutive sites is governed by a “soft” simple exclusion principle and by attracting or repelling forces between neighboring particles. Using tools from contraction theory, we prove that the model admits a unique steady-state and that every trajectory converges to this steady-state. Analysis and simulations of the effect of the attracting and repelling forces on this steady-state highlight the crucial role that these forces may play in increasing the steady-state flow, and reveal that this increase stems from the alleviation of traffic jams along the lattice. Our theoretical analysis clarifies microscopic aspects of complex multi-particle dynamic processes.

## Introduction

Biological processes are governed by complex interactions between multiple particles that are confined in special compartments [1]. One of the most important examples of such processes is biological intracellular transport, which is carried by motor proteins (e.g., kinesins, dyneins, and myosins) [2]. These motor proteins, which are also known as biological molecular motors, can catalyze the reaction of adenosine triphosphate (ATP) hydrolysis, while at the same time converting the energy produced during this chemical reaction into a mechanical work required for their movements along cellular filaments (such as microtubules and actin filaments) [2].

Experimental observations clearly show that motor proteins usually function in large groups, suggesting that the interactions between the motors cannot be ignored [3, 4]. Understanding the collective behavior of molecular motors is critical for uncovering mechanisms of complex biological processes [2, 5, 6]. From a theoretical point of view, intracellular transport processes are usually described using non-equilibrium multi-particle lattice models [3]. In these models, the molecular motors are typically represented by particles that hop along the lattice, and the lattice sites model the binding locations of the motors along the filaments (or tracks). For a general review on transport and traffic phenomena in biological systems see for example [2, 5, 3, 4].

A standard model from non-equilibrium statistical mechanics for molecular motors traffic (and numerous other processes) is the totally asymmetric simple exclusion process (TASEP) [7, 8, 9]. This is also the standard model for ribosome flow during mRNA translation (see, e.g. [10, 11, 8]). In TASEP, particles hop randomly along a unidirectionally ordered lattice of sites. Simple exclusion means that a particle cannot move into a site that is already occupied by another particle, and thus each site can be either empty or occupied by a single particle. This models moving biological particles like ribosomes and motor proteins that have volume and thus cannot overtake a moving particle in front of them. This hard exclusion principle creates an intricate indirect coupling between the particles. In particular, a slowly moving particle may lead to the formation of a traffic jam behind it.

To describe moving biological molecules with large sizes, a version of TASEP with extended objects has been introduced and analyzed [12, 13, 14]. In this model, each particle covers  lattice sites. Thus a particle occupies sites  for some , and it can hop to site  provided that site  is empty. This is used, for example, for modeling mRNA translation as it is known that every ribosome (the particle) covers several codons (sites) along the mRNA molecule [13].

There exist two versions of TASEP that differ by their boundary conditions. In TASEP with open boundary conditions the two sides of the chain are connected to two particle reservoirs with constant concentrations, and the particles can hop into the lattice chain (if the first site is empty) and out of the chain (if the last site is occupied). In the open boundary homogeneous TASEP (HTASEP), all the transition rates within the lattice are assumed to be equal and normalized to one, and thus the model is specified by an input rate , an exit rate , and a parameter denoting the number of sites along the lattice. In TASEP with periodic boundary conditions the chain is closed into a ring, and a particle that hops from the last site returns to the first site. TASEP has been widely utilized for studying various natural and artificial processes, including vehicular traffic flow, mRNA translation, surface growth, communication networks, and more [3, 15].

Ref. [16] used HTASEP with periodic boundary conditions to analyze transport on a lattice in the presence of local interactions between particles and substrate, illustrating the effect of local conformation of the substrate on the characteristics of the flow of molecular motors. TASEP with particle interactions and with periodic boundary conditions was studied in [17], and with open boundary conditions in [18, 19, 20, 21]. Specifically, the authors in [20, 21] proposed a modified TASEP model that incorporates the realistic observed feature of nearest-neighbor interactions. In this model, the transition rate in every site along the lattice depends on the states of four consecutive sites. Their conclusions were that weak repulsive interaction results in maximal flux, and that the molecular motors are influenced more strongly by attractive interactions.

Unfortunately, rigorous analysis of TASEP is non-trivial, and exact solutions exist only in special cases, for example when considering the model with the homogeneous rates (HTASEP). Typically, the non-homogeneous case and cases that include other local interactions are only studied via various approximations and extensive Monte Carlo computer simulations. These simulations are run until convergence to a (stochastic) steady-state, yet without a rigorous proof that convergence indeed takes place for all the feasible parameter values.

In this paper, we introduce a new deterministic model for the flow of motor proteins along a one-dimensional lattice of sites with nearest-neighbor interactions between the motors. The flow of the motor proteins is unidirectional, and it satisfies a “soft” simple exclusion principle. The nearest-neighbor effect is modeled by two “force” interactions with parameters  and . It is more convenient to explain the effect of these interactions in “particle-like” terms, although in the new model the density in every site takes values in the range  (and not ).

Consider a transition of a particle from site  to site . If site  is already occupied then the rate of movement depends on a parameter  that represents an “attachment/detachment force” when generating new neighbors. A value  [] means that the particle will tend [not] to hop forward, as there is a strong attraction [repulsion] to the particle in site . On the other-hand, if site  is already occupied then the rate of movement depends on a parameter  that represents an “attachment/detachment force” when breaking from old neighbors. A value  [] means that the particle will tend [not] to hop forward, as there is a strong repulsion [attraction] from the neighboring particle in site . A value of  [] implies no attachment/detachment force when generating new neighbors [when breaking from old neighbors].

An important advantage of our model is that it is highly amenable to rigorous analysis even for non-homogenous transition rates. We prove, for example, that the dynamics always converges to a steady-state density along the lattice. Thus, the flow also converges to a steady-state value. This steady-state depends on the lattice size, the transition rates, and the parameters , , but not on the initial density along the lattice (i.e. the initial conditions). Analysis and simulations of the effect of the attracting and repelling forces on this steady-state highlight the crucial role that these forces may play in increasing the steady-state flow, and reveal that this increase stems from the alleviation of traffic jams along the lattice. It is well-known that molecular motors indeed form traffic jams and that these have important biological implications (see, e.g. [22, 23, 24]). In particular, analysis and simulations of the model reveal a new regime that may be interpreted as the “opposite” of a traffic jam along the lattice.

Our approach extends a deterministic mathematical model that has been used for describing and analyzing the flow of ribosomes along the mRNA molecule during the process of mRNA translation. The next section provides a brief overview of this model.

### The Ribosome Flow Model (RFM)

The RFM [25] is a nonlinear, continuous-time, compartmental model for the unidirectional flow of “material” along a one-dimensional chain of consecutive compartments. It can be derived via a mean-field approximation of TASEP with open boundary conditions [3, Section 4.9.7] [7, p. R345]. The RFM includes parameters: controls the initiation rate, the exit rate, and , , the transition rate from site  to site . The state variable , , describes the normalized amount of “material” (or density) at site  at time , where  [] indicates that site  is completely full [completely empty] at time . Thus, the vector  describes the density profile along the chain at time . The output rate at time is  (see Fig. 1).

Let , and . The dynamics of the RFM with sites is given by the following set of nonlinear ODEs:

 ˙xi=λi−1xi−1(1−xi)−λixi(1−xi+1),i=1,…,n. (1)

This can be explained as follows. The flow of material from site  to site  at time  is . This flow increases with the density at site , and decreases as site  becomes fuller. This corresponds to a “soft” version of a simple exclusion principle. Note that the maximal possible flow from site  to site  is the transition rate . Thus Eq. (1) simply states that the change in the density at site at time is the input rate to site (from site ) at time minus the output rate (to site ) at time .

The trajectories of the RFM evolve on the compact and convex state-space

 Cn:={x∈Rn:xi∈[0,1],i=1,…,n}.

Let [] denote the interior [boundary] of . Ref. [26] has shown that the RFM is a tridiagonal cooperative dynamical system [27], and consequently (1) admits a unique steady-state density  that is globally asymptotically stable, that is, for all (see also [28]). This means that trajectories corresponding to different initial conditions all converge to the same steady-state density . In particular, the density at the last site  converges to the value , so the output rate  converges to a steady-state value .

An important advantage of the RFM (e.g. as compared to TASEP) is that it is amenable to mathematical analysis using tools from systems and control theory. Furthermore, most of the analysis hold for the general, non-homogeneous case (i.e. the case where the transition rates  differ from one another). For more on the analysis of the RFM and its biological implications, see [29, 30, 26, 31, 28, 32, 33, 34, 35, 36, 37].

In this paper, we extend the RFM to include nearest-neighbor interactions, namely, binding and repelling actions that are dynamically activated for each site based on the state of its neighboring sites. A parameter [] controls the binding/repelling forces between two existing [new] neighbors. We refer to the new model as the excluded flow with local repelling and binding model (EFRBM). It is important to note that this is significantly different from the RFM. For example, the EFRBM, unlike the RFM, is not a cooperative system [27]. Also, in the RFM the dynamics at site is directly affected by its two nearest neighbors sites, whereas in the EFRBM the dynamics is directly affected by the density in four neighboring sites. Thus, unlike the RFM, the EFRBM is not a tridiagonal system. Also, the RFM has been used to model ribosome flow, whereas here we apply the EFRBM to study the flow of motor proteins.

We show that the EFRBM is a contractive dynamical system. This holds for any set of feasible transition rates and local interaction forces including the case of non-homogeneous transition rates. This implies that the EFRBM admits a unique steady-state that is globally asymptotically stable. Thus, every set of parameters corresponds to a unique steady-state output rate. We analyze the behavior of this steady-state under the assumption  that follows from fundamental thermodynamic arguments (see [38]). We show that a small neighbor-repelling force (i.e. small  and thus a large ) leads to a small output rate. Analysis and simulations show that this is due to the formation of traffic jams at the beginning of the lattice. On the other-hand, a strong neighbor-repelling force (i.e. large  and small ) lead to a high output rate. In this case, an interesting phenomena emerges: the density in every second site goes to zero. This “separation of densities” is the “opposite” of a traffic jam. These results highlight the impact of traffic jams on the output rate.

The remainder of this paper is organized as follows. The next section describes the EFRBM. The following two sections describe our main analysis results and their biological implications. This includes analysis of the asymptotic behavior of the EFRBM, and the effects of the nearest-neighbor interactions on the steady-state behavior of the EFRBM. The final section summarizes and describes several directions for further research. To increase the readability of this paper, all the proofs are placed in the Appendix.

## The EFRBM

The EFRBM with sites includes  parameters:

• , , controls the transition rate from site to site , where [] controls the input [output] rate.

• is the attachment/detachment force between any two existing (consecutive) neighbors.

• is the attachment/detachment force between any two new (consecutive) neighbors.

Fig. 2 depicts the four possible transition scenarios from site  to site , and the rates in each case. For simplicity, we use a schematic “particle-like” explanation, although in the EFRBM the state-variables represent a normalized material density in the range  and not a binary choice  like in TASEP. If both sites and do not contain particles, the transition rate is simply , as in the RFM. If a particle is located at site [] but site  [] is empty then the transition rate is  []. If both sites contain particles the transition rate is .

The EFRBM also includes  state-variables , . Just like in the RFM, describes the normalized density at site  at time , where  [] means that the site is completely empty [full].

To state the dynamical equations describing the EFRBM we introduce more notation. Let , , and denote

 zi(t):={xi(t),i=1,…,n,0,otherwise. (2)

Then the EFRBM is described by

 ˙xi=gi−1(x)−gi(x),i=1,…,n, (3)

where

 gi(x):=λixi(1−xi+1)(1+(q−1)zi+2)(1+(r−1)zi−1). (4)

We now explain these equations. The term  represents the flow from site  to site , so Eq. (3) means that the change in the density at site  is the inflow from site minus the outflow to site . To explain Eq. (4), consider for example the case  (and assume that ). Then (4) yields

 g2(x)=λ2x2(1−x3)(1+(q−1)x4)(1+(r−1)x1). (5)

The term  means that the flow from site  to site  increases with the density at site . The term  represents soft exclusion: as the density at site  increases, the transition from site  to site  gradually decreases. The term  represents the fact that the flow into site  also depends on the density at site : if  [] then the transition increases [decreases] with , that is, the “particles” at site  “attract” [“repel”] the particles that move from site  to site . The term  is similar but represents an attachment/detachment force between the “particles” in sites  and .

Note that for , , and thus in this case the EFRBM reduces to the RFM (see (1)). On the other hand, if then . This represents a kind of an “extended objects” RFM, as the transition from site  to site  decreases with the density in sites , , and .

###### Remark 1

It is useful to think of the EFRBM as an RFM with time-varying transition rates. For example, we can write (5) as

 g2(x(t))=η2(t)x2(t)(1−x3(t)),

where . Note that this time-varying transition rate depends on  (i.e., the fixed site to site transition rate), and also on and and the time-varying densities in the neighboring sites, as these determine the interaction forces between the moving particles.

We denote the flow from site  to the environment by

 R(t) :=λnxn(t)(1+(r−1)xn−1(t)). (6)

This is the output rate at time .

###### Example 1

The EFRBM with  sites is given by:

 ˙x1 =λ0(1−x1)(1+(q−1)x2)−λ1x1(1−x2)(1+(q−1)x3), ˙x2 =λ1x1(1−x2)(1+(q−1)x3)−λ2x2(1−x3)(1+(r−1)x1), ˙x3 =λ2x2(1−x3)(1+(r−1)x1)−λ3x3(1+(r−1)x2). (7)

If  then this becomes

 ˙x1 =λ0(1−x1)(1−x2)−λ1x1(1−x2)(1−x3), ˙x2 =λ1x1(1−x2)(1−x3)−λ2(1−x1)x2(1−x3), ˙x3 =λ2(1−x1)x2(1−x3)−λ3(1−x2)x3. (8)

On the other-hand, for and  (1) becomes

 ˙x1 =λ0(1−x1)−λ1x1(1−x2), ˙x2 =λ1x1(1−x2)−λ2(1−x1)x2(1−x3), ˙x3 =λ2(1−x1)x2(1−x3)−λ3(1−x2)x3, (9)

Following [38] (see also [39]), we view creating and breaking a pair of particles as opposite chemical transitions, so by detailed balance arguments: , where  is the interaction energy. As in [38], we also assume that  is equally split between the creation and breaking processes, so

 q=exp(E2KBT),r=exp(−E2KBT). (10)

This has a clear physical meaning. If  the interaction is attractive, so the particle moves faster when creating a new pair  since the energy of the system decreases by . On the other-hand, breaking out of the cluster increases the energy by  and the transition rate is thus slowed down (). Similarly, the case  corresponds to a repulsive interaction and then  and . Note that (10) implies in particular that

 rq=1. (11)

In this case, the EFRBM contains parameters: , and (as ). Note that if (11) holds then (4) becomes

 gi(x)=λixi(1−xi+1)(1−r−1rzi+2)(1+(r−1)zi−1). (12)

The next section derives several theoretical results on the dynamical behavior of the EFRBM. Recall that all the proofs are placed in the Appendix.

## Asymptotic behavior of the EFRBM

Let  denote the solution of the EFRBM at time  for the initial condition .

### Invariance and persistence

The next result shows that the -dimensional unit cube  is an invariant set of the EFRBM, that is, any trajectory that emanates from an initial condition in  remains in  for all time. Furthermore, any trajectory emanating from the boundary of  “immediately enters” . This is a technical result, but it is important as in the interior of  the EFRBM admits several useful properties.

###### Proposition 1

Assume that . For any  there exists such that

 d≤xi(t+τ,a)≤1−d,

for all , all , and all .

This means that all the trajectories of the EFRBM enter and remain in the interior of  after an arbitrarily short time. In particular, both and are invariant sets of the EFRBM dynamics.

From a biological point of view this means that if the system is initiated such that every density is in  then this remains true for all time , so the equations “make sense” in this respect. Furthermore, after an arbitrarily short time the densities are all in , i.e. any completely empty [full] site immediately becomes not completely empty [full].

### Contraction

Differential analysis and in particular contraction theory proved to be a powerful tool for analyzing the asymptotic behavior of nonlinear dynamical systems. In a contractive system, trajectories that emanate from different initial conditions approach each other at an exponential rate [40, 41, 42].

For our purposes, we require a generalization of contraction with respect to (w.r.t.) a fixed norm that has been introduced in [43]. Consider the time-varying dynamical system:

 ˙x(t)=f(t,x(t)), (13)

whose trajectories evolve on an invariant set  that is compact and convex. Let  denote the solution of (13) at time  for the initial condition . The dynamical system (13) is said to be contractive after a small overshoot (SO) [43] on  w.r.t. a norm  if for any  there exists  such that

 |x(t,t0,a)−x(t,t0,b)|≤(1+ε)exp(−(t−t0)ℓ)|a−b|,

for all and all . Intuitively speaking, this means that any two trajectories of the system approach each other at an exponential rate , but with an arbitrarily small overshoot of .

Let  denote the  norm, i.e. for .

###### Proposition 2

The EFRBM with is SO on  w.r.t. the  norm, that is, for any  there exists  such that

 |x(t,a)−x(t,b)|1≤(1+ε)exp(−ℓt)|a−b|1, (14)

for all and all .

From a biological point of view this means the following. The state of the system at any time  is a vector describing the density at each site at time . We measure the distance between any two density vectors using the  vector norm. Suppose that we initiate the system with two different densities. This generates two different solutions of the dynamical system. The distance between these solutions decreases with time at an exponential rate.

The next example demonstrates this contraction property. Let  [] denote the column vector of  ones [zeros].

###### Example 2

Consider the EFRBM with dimension , and parameters , , , , and . Fig. 3 depicts , with and , as a function of time for . It may be seen that the distance between the two trajectories goes to zero at an exponential rate.

Prop. 2 implies that the EFRBM satisfies several important asymptotic properties. These are described in the following subsections.

### Global asymptotic stability

Write the EFRBM (3) as . Since the compact and convex set  is an invariant set of the dynamics, it contains at least one steady-state. That is, there exists  such that . By Proposition 1. Using (14) with yields the following result.

###### Corollary 1

Assume that . Then the EFRBM admits a unique steady-state  that is globally asymptotically stable, i.e.

 limt→∞x(t,a)=e, for all a∈Cn.

This means that any solution of the EFRBM converges to a unique steady-state density (and thus a unique steady-state output rate) that depends on the rates , and the parameters  and , but not on the initial condition. From a biological point of view, this means that the system always converges to a steady-state density and a corresponding steady-state output rate, and thus it makes sense to study how these depend on the various parameters.

Note that the assumption that  cannot be dropped. Indeed, Eq. (1), corresponding to a EFRBM with , and , admits a continuum of steady-states.

###### Example 3

Fig. 4 depicts the trajectories of (3) with , , , , , and , for several initial conditions. It may be seen that all trajectories converge to a unique steady-state . (All the numerical values in the simulations described in this paper are to four digit accuracy.)

The rigorous proof that every trajectory converges to a steady-state is important, as it implies that after some time the densities are very close to their steady-state values. The next step is to analyze this steady-state density and the corresponding steady-state output rate, and explore how these are related to the various parameters of the model.

At steady-state, (i.e. for ) the left-hand side of all the equations in (3) is zero (i.e. , ), so for all . This implies that

 λ0 (1−e1)(1+(q−1)e2) =λ1e1(1−e2)(1+(q−1)e3) =λ2e2(1−e3)(1+(q−1)e4)(1+(r−1)e1) =λ3e3(1−e4)(1+(q−1)e5)(1+(r−1)e2) ⋮ (15) =λn−1en−1(1−en)(1+(r−1)en−2) =λnen(1+(r−1)en−1),

and also that the steady-state flow satisfies

 R=gi(e),i=0,…,n. (16)

In particular,  and since every , the steady-state flow is positive (i.e. a left-to-right flow) for any .

Also, for the case  it follows from  that for , , whereas for  it follows from  that , so

 R≤max{λ0,λn}.

This means in particular that the output rate is always bounded.

###### Fact 1

It follows from (Analysis of the steady-state) that if we multiply all the s by a parameter then will not change, i.e. . Thus, by (16) , for all , that is, the steady-state flow [density] is homogeneous of degree one [zero] w.r.t. the s.

In the spacial case  the steady-state equations (Analysis of the steady-state) can be solved in closed-form.

###### Fact 2

Consider the EFRBM with  and . Define

 a1 :=(1−1r)(λ2r+λ1)+λ1λ2λ0. (17)

Then  is given by

 e2 =λ1+λ2+a1−√(λ1+λ2+a1)2−4a1λ12a1, e1 =λ2e2λ1+(λ2(1−r)−λ1)e2. (18)

Note that even in this case the expression for  is non-trivial.

Let  denote the set of  dimensional vectors with all entries positive. Let  denote the set of parameters in the EFRBM with dimension . The results above imply that there exists a function  such that  is the unique steady-state of the EFRBM with parameters .

###### Proposition 3

The function  is analytic.

This result allows in particular to consider the derivatives of the steady-state density  and the steady-state output rate  w.r.t. small changes in some of the parameters , that is, the sensitivity of the steady-state w.r.t. small changes in the parameters.

## Effect of nearest-neighbor interactions

We begin with several simulations demonstrating the effect of the parameter  (and ) on the steady-state of the EFRBM.

###### Example 4

Consider a EFRBM with  and rates . Fig. 5 depicts the steady-state output rate  as a function of . It may be seen that  monotonically increases with . In particular, for  (i.e., the RFM) , wheres for , , that is, the steady-state flow is increased by about . When considering the comparison with the RFM, one should bear in mind that the EFRBM corresponds to an RFM with time-varying rates  that may effectively be much higher than the fixed rates . We assume that the energy that is needed to generate these higher rates comes from the additional interaction forces between the particles.

The next example demonstrates that the increase in  as  increases is because the neighbor-repelling forces lead to an alleviation of traffic jams.

###### Example 5

Consider the EFRBM with dimension , , , , , , , and . Consider first the case (i.e., the RFM). The steady-state density is:

 e=[0.84430.84630.79560.95100.18140.1416]′,

and the corresponding steady-state flow is . Note that since  is high and  is low, , indicating a traffic jam at site . Consider now the case (i.e. ). The steady-state density is now

and the corresponding steady-state output is . Note that now the density at site decreased relative to the case, and that . Note also that . This means that the introduction of a “neighbor-repelling” force (i.e. ) alleviated the traffic jam, reduced the total steady-state occupancy, and increased the steady-state flow.

Fig. 6 depicts the steady-state densities in this example as a function of . It may be observed that , , monotonically decreases with , and that slightly increases with . Note that since the occupancy at site is not affected by , but only by , increasing should indeed increase

### Extreme interactions

To gain more insight on the effect of the nearest-neighbor interactions on the steady-state behavior, it is useful to consider the cases when  (so ) and  (so ).

#### The case r→0

Intuitively speaking, a low value of  corresponds to: (1) a strong attachment between existing nearest neighbors (small ); and (2) a high tendency for moving forward if this involves creating new neighbors (large ). As we will see this leads to the formation of traffic jams and, consequently, to a sharp decrease in the output rate.

###### Example 6

Consider a EFRBM with dimension  and rates , . For (recall that ), the steady-state values are:

 e=[0.99080.98990.90620.89780.98410.5678]′,R=0.0913.

For , the steady-state values are:

 e=[0.99980.99990.99010.99000.98990.4970]′,R=0.0099.

For , the steady-state density values are:

 e=[0.99960.99990.99500.99500.99500.4986]′,R=0.0050.

Fig. 7 depicts the steady-state values for the three  values. It may be observed that as  decreases the density in the first five sites increases to one, i.e. these sites become completely full, and the output rate goes to zero. Note that this highlights the negative effect of traffic jams on the output rate.

We now rigorously analyze the case  for the EFRBM with  and .

###### Example 7

Consider the EFRBM with  and . Expanding  and  in (2) as a Taylor series in  yields

 e2=1−λ2λ1r+o(r),e1=1+o(r), (19)

where every  denotes a function  satisfying . Thus, . This implies in particular that

 limr→0e1=limr→0e2=1,limr→0R=0.

Thus, when , both steady-state densities go to one,111Note that (19) implies that  goes to one faster than . that is, the sites become completely full, and consequently the steady-state output rate goes to zero.

The next result analyzes the case .

###### Proposition 4

The steady-state densities in the EFRBM with  satisfy

 e1(r) =1−λ2λ3λ0(λ2+λ3)r2+o(r2), e2(r) =1−λ3λ1r2+o(r2), e3(r) =λ2λ2+λ3+o(r), (20)

and

 R(r)=λ2λ3λ2+λ3r+o(r). (21)

Note that this implies that

 limr→0e1(r)=limr→0e2(r)=1, and limr→0R(r)=0,

so again as  sites at the beginning of the lattice become completely full and consequently the output rate goes to zero.

Summarizing, as goes to the repelling force between existing neighbors is very weak, and the binding force when forming new neighbors is very strong, leading to the formation of traffic jams at the beginning of the lattice. Consequently, the steady-state flow goes to zero.

We now turn to consider the opposite case, that is, .

#### The case r→∞

A large value of  corresponds to: (1) strong repulsion between existing nearest neighbors (large ); and (2) a low tendency for moving forward if this involves creating new neighbors (small ). As we will see below, this leads to a phenomena that may be regarded as the opposite of traffic jams, that is, a complete “separation of the densities” along the lattice.

###### Example 8

Consider the EFRBM with sites and rates For (recall that ),

 e=[0.73220.69500.51830.45580.46570.2329]′,R=0.2678.

For ,

 e=[0.52620.00150.24980.00220.20070.0020]′,R=0.4729.

For ,

 e=[0.52610.00010.24950.00020.20010.0002]′,R=0.4734.

Fig. 8 depicts these steady-state values for the three values. Note that the steady-state values for  and  cannot be distinguished. It may be observed that the values , , decrease to zero as  increases. In other words, in every pair of consecutive sites one density is very small. This “separation of densities” represents the opposite of a traffic jam. This leads to a substantial increase in the output rate  as  increases.

We now rigorously analyze the case  for the EFRBM with  and .

###### Example 9

Consider the EFRBM with . Expanding  in (2) as a Taylor series in  yields

 e2=λ1λ2q+o(q),e1=λ0λ0+λ1+o(q),

so

 limr→∞e2=0,limr→∞e1=λ0λ0+λ1,limr→∞R=λ0λ1λ0+λ1.

Thus, in this case the density at site  goes to zero, and this yields a positive steady-state output rate.

###### Proposition 5

The steady-state densities in the EFRBM with  satisfy

 e1(q) =a1+b1q+o(q), e2(q) e3(q) =a3+b3q+o(q), (22)

with , and

 R(q)=λ0(1−a1)+λ0((a1−1)λ1λ2−b1)q+o(q). (23)

Note that this implies that

 limr→∞e2(r)=0, and limr→∞R(r)>0,

so again as  the density at site  goes to zero and the output rate is positive.

## Discussion

Motor proteins and other moving biological particles interact with their neighbors. Indeed, it is known that cellular cargoes are often moved by groups of motor proteins, and recent findings suggest that the bounding time of kinesins on microtubules depend on the presence of neighbors.

To study the effect of such interactions, we introduced a new deterministic compartmental model, the EFRBM, for the flow of particles along an ordered lattice of sites where the transition rates between sites depend both on properties of the lattice and on nearest-neighbor interactions between the particles. The properties of the lattice are modeled using transition rates  between sites. The nearest-neighbor interactions between the particles are modeled using two parameters:  that represents the tendency of a moving particle to break from an existing neighbor, and  that represents the tendency of a particle to move into a site such that it forms new neighbors (see Fig. 2).

The EFRBM is based on a mean-field ansatz neglecting high-order correlations of occupations between neighboring sites. It is possible to use our framework also to derive a more complete model based on binary occupation densities and transitions described by a continuous-time master equation (see, e.g. the interesting paper [44] in which this was done for granular channel transport). However, in such a model the state-variables at time  represent the probability of each configuration at time , and the number of possible configurations grows exponentially with the number of sites . On the other-hand, the EFRBM includes  (nonlinear) ODEs for  sites. Another important advantage of the EFRBM is that it is amenable to analysis using tools from systems and control theory, even in the non-homogeneous case. This allows to rigorously study, for example, the effect of the nearest-neighbors interactions on the steady-state behavior of the EFRBM for any set of transition rates. Our results show that suitable forces between nearby particles can greatly increase the output rate, and reveal that the underlying mechanism for this is the alleviation of traffic jams along the lattice. In particular, when the parameter  is very large and  is very small, the steady-state density is such that any second site is empty. This represents the “opposite” of a traffic jam, and increases the steady-state flow.

The phenomenological model introduced here may prove useful for other applications as well. For example, an important problem in vehicular traffic is to understand how human drivers react to nearby cars. One may also consider implementing appropriate nearest-neighbor dynamics in algorithms that control autonomous vehicles in order to reduce traffic jams and increase the flow. Of course, implementing this with a very large  (or ) means very high effective transition rates, but our results suggest that even for  not much larger than one the increase in the flow is non-negligible. Another interesting topic for further research is generalizing the EFRBM to include the possibility of attachment/detachment of particles from intermediate sites in the lattice (see [45] for some related ideas).

## Acknowledgments

We thank the anonymous referees and the Associate Editor for their helpful comments.

The research of YZ is partially supported by the Edmond J. Safra Center for Bioinformatics at Tel Aviv University. The research of MM is partially supported by research grants from the Israeli Ministry of Science, Technology & Space, the US-Israel Binational Science Foundation, and the Israeli Science Foundation.

## Appendix: Proofs

Proof of Proposition 1. The fact that  is an invariant set of the dynamics follows immediately from the equations of the EFRBM. Let

 ηi(t):=λi(1+(q−1)zi+2(t))(1+(r−1)zi−1(t)),i=0,…,n, (24)

with the s defined in (2). By (3), the EFRBM can be written as

 ˙xi(t)=ηi−1(t)xi−1(t)(1−xi(t))−ηi(t)xi(t)(1−xi+1(t)). (25)

This is just the RFM (see (1)), but with time-varying rates . Let , and . It follows from (24) that  for all  and for all . Note that for  every  is strictly positive. In other words, all the time-varying rates are uniformly separated from zero and uniformly bounded. Now the proof of Proposition 1 follows from the results in [28]

Proof of Proposition 2. Combining the representation in (25) with the uniform boundedness of the rates, Proposition 1, and the results in [43] imply that the EFRBM is contractive after a small overshoot and short transient (SOST) on . Also, Proposition 4 in [43] implies that for the EFRBM the properties of SOST and SO are equivalent, and this completes the proof.

Proof of Fact 2. Consider the EFRBM with  and . Then (Analysis of the steady-state) becomes

 λ0(1−e1)(1+(1r−1)e2) =λ1e1(1−e2) =λ2e2(1+(r−1)e1).

This yields

 e1=λ2e2λ1+(λ2(1−r)−λ1)e2 (26)

and

 a1e22+a2e2+λ1=0,

with  defined in (17) and . The feasible solution (i.e. the one satisfying  for any set of parameter values) is given by

 e2=−a2−√a22−4a1λ12a1,

and (26).

Proof of Prop. 3. To emphasize the dependence on the parameters, write the EFRBM as , where . Note that  is an analytic function. Then the steady-state satisfies the relation . The Jacobian matrix of this relation with respect to  is

 J(x;v):=∂∂xf(x;v),

which is just the Jacobian of the dynamics. Fix  and let  denote the corresponding steady-state, that is,  and . Suppose that there exists a matrix measure  such that . This implies in particular that  is Hurwitz (see e.g. [46]), so it is not singular and invoking the implicit function theorem implies that the mapping  is analytic. It follows from the results in [28] that such a matrix measure  indeed exists, and this completes the proof.

Proof of Prop. 4. Expand , , as

 ei=ai+bir+cir2+o(r2). (27)

Recall that the steady-state equations are given by , with the s given in (4). Substituting (27) yields

 g0(e) =λ0(1−a1)a2r+…, g1(e) =λ1(1−a2)a3r+…, g2(e) =λ2(a1−1)a2(a3−1)+…, g3(e) =λ3(1−a2)a3+….

Since  is bounded, we conclude that

 (1−a1)a2=(1−a2)a3=0. (28)

Assume for the moment that . Then

 g2(e) =λ2(a1−1)(a3−1)b2r+…, (29) g3(e) =λ3a3+…,

and this implies that . Now,