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

TL;DR

This paper introduces a deterministic model for one-dimensional particle flow with local interactions, revealing how attractive and repulsive forces influence steady-state flow and traffic congestion.

Contribution

The paper develops a novel deterministic compartmental model incorporating local interactions and proves convergence to a unique steady-state using contraction theory.

Findings

Attracting forces increase steady-state flow.

Repelling forces help alleviate traffic jams.

Model predictions align with simulation results.

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.

Equations

$${\dot{x}}_i={\lambda }_{i-1}{x}_{i-1}(1-x_i)-{\lambda }_i{x}_i(1-x_{i+1}),i=1,\dots ,n.$$

$$C^n:=\{x\in {\mathbb{R}}^n:x_i\in [0,1],i=1,\dots ,n\}.$$

$$z_i(t):=\{\begin{array}{ll}{x}_i(t),& i=1,\dots ,n,\\ 0,& \text{otherwise}.\end{array}$$

$${\dot{x}}_i=g_{i-1}(x)-g_i(x),i=1,\dots ,n,$$

$$g_i(x):={\lambda }_i{x}_i(1-x_{i+1})(1+(q-1)z_{i+2})(1+(r-1)z_{i-1}).$$

$$g_2(x)={\lambda }_2{x}_2(1-x_3)(1+(q-1)x_4)(1+(r-1)x_1).$$

$$g_2(x(t))={\eta }_2(t)x_2(t)(1-x_3(t)),$$

$$R(t)$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{x}}_3$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{x}}_3$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{x}}_3$$

$$q=\exp \left(\frac{E}{2K_{B}T}\right),r=\exp \left(\frac{-E}{2K_{B}T}\right).$$

$$rq=1.$$

$$g_i(x)={\lambda }_i{x}_i(1-x_{i+1})(1-\frac{r-1}{r}{z}_{i+2})(1+(r-1)z_{i-1}).$$

$$d\le x_i(t+\tau ,a)\le 1-d,$$

$$\dot{x}(t)=f(t,x(t)),$$

$$|x(t,t_0,a)-x(t,t_0,b)|\le (1+\epsilon )\exp (-(t-t_0)\ell )|a-b|,$$

$$|x(t,a)-x(t,b){|}_1\le (1+\epsilon )\exp (-\ell t)|a-b{|}_1,$$

$$\underset{t\to \infty }{\lim }x(t,a)=e,\text{ for all }a\in C^n.$$

$${\lambda }_0$$

$$={\lambda }_1{e}_1(1-e_2)(1+(q-1)e_3)$$

$$={\lambda }_2{e}_2(1-e_3)(1+(q-1)e_4)(1+(r-1)e_1)$$

$$={\lambda }_3{e}_3(1-e_4)(1+(q-1)e_5)(1+(r-1)e_2)$$

$$⋮$$

$$={\lambda }_{n-1}{e}_{n-1}(1-e_n)(1+(r-1)e_{n-2})$$

$$={\lambda }_n{e}_n(1+(r-1)e_{n-1}),$$

$$R=g_i(e),i=0,\dots ,n.$$

$$R\le \max \{{\lambda }_0,{\lambda }_n\}.$$

$$a_1$$

$$e_2$$

Full text

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

Yoram Zarai1, Michael Margaliot2,∗, Anatoly B. Kolomeisky3

**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.

$\ast$ 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 $\ell>1$ lattice sites. Thus a particle occupies sites $i,\dots,i+\ell-1$ for some $i$, and it can hop to site $i+1$ provided that site $i+\ell$ 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 $\alpha$, an exit rate $\beta$, and a parameter $N$ 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 $q$ and $r$. 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 $[0,1]$ (and not $\{0,1\}$).

Consider a transition of a particle from site $i$ to site $i+1$. If site $i+2$ is already occupied then the rate of movement depends on a parameter $q\geq 0$ that represents an “attachment/detachment force” when generating new neighbors. A value $q>1$ [$q<1$] means that the particle will tend [not] to hop forward, as there is a strong attraction [repulsion] to the particle in site $i+2$. On the other-hand, if site $i-1$ is already occupied then the rate of movement depends on a parameter $r\geq 0$ that represents an “attachment/detachment force” when breaking from old neighbors. A value $r>1$ [$r<1$] means that the particle will tend [not] to hop forward, as there is a strong repulsion [attraction] from the neighboring particle in site $i-1$. A value of $q=1$ [$r=1$] 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 $q$, $r$, 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 $n$ 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 $n+1$ parameters: $\lambda_{0}>0$ controls the initiation rate, $\lambda_{n}>0$ the exit rate, and $\lambda_{i}>0$, $i=1,\dots,n-1$, the transition rate from site $i$ to site $i+1$. The state variable $x_{i}(t):\mathbb{R}_{+}\to[0,1]$, $i=1,\dots,n$, describes the normalized amount of “material” (or density) at site $i$ at time $t$, where $x_{i}(t)=1$ [$x_{i}(t)=0$] indicates that site $i$ is completely full [completely empty] at time $t$. Thus, the vector $x(t):=\begin{bmatrix}x_{1}(t)&\dots&x_{n}(t)\end{bmatrix}^{\prime}$ describes the density profile along the chain at time $t$. The output rate at time $t$ is $R(t):=\lambda_{n}x_{n}(t)$ (see Fig. 1).

Let $x_{0}(t)\equiv 1$, and $x_{n+1}(t)\equiv 0$. The dynamics of the RFM with $n$ sites is given by the following set of $n$ nonlinear ODEs:

[TABLE]

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

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

[TABLE]

Let $\operatorname{{\mathrm{I}nt}}(C^{n})$ [$\partial C^{n}$] denote the interior [boundary] of $C^{n}$. Ref. [26] has shown that the RFM is a tridiagonal cooperative dynamical system [27], and consequently (1) admits a unique steady-state density $e=e(\lambda_{0},\dots,\lambda_{n})\in\operatorname{{\mathrm{I}nt}}(C^{n})$ that is globally asymptotically stable, that is, $\lim_{t\to\infty}x(t,a)=e$ for all $a\in C^{n}$ (see also [28]). This means that trajectories corresponding to different initial conditions all converge to the same steady-state density $e$. In particular, the density at the last site $x_{n}(t)$ converges to the value $e_{n}$, so the output rate $R(t)$ converges to a steady-state value $R:=\lambda_{n}{e}_{n}$.

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 $\lambda_{i}$ 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 $r$ [$q$] 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 $i$ 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 $rq=1$ that follows from fundamental thermodynamic arguments (see [38]). We show that a small neighbor-repelling force (i.e. small $r$ and thus a large $q=1/r$) 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 $r$ and small $q$) 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 $n$ sites includes $n+3$ parameters:

• •

$\lambda_{i}>0$, $i=0,\dots,n$, controls the transition rate from site $i$ to site $i+1$, where $\lambda_{0}$ [$\lambda_{n}$] controls the input [output] rate.

• •

$r\geq 0$ is the attachment/detachment force between any two existing (consecutive) neighbors.

• •

$q\geq 0$ is the attachment/detachment force between any two new (consecutive) neighbors.

Fig. 2 depicts the four possible transition scenarios from site $i$ to site $i+1$, 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 $[0,1]$ and not a binary choice $\{0,1\}$ like in TASEP. If both sites $i-1$ and $i+2$ do not contain particles, the transition rate is simply $\lambda_{i}$, as in the RFM. If a particle is located at site $i-1$ [$i+2$] but site $i+2$ [$i-1$] is empty then the transition rate is $\lambda_{i}r$ [$\lambda_{i}q$]. If both sites contain particles the transition rate is $\lambda_{i}rq$.

The EFRBM also includes $n$ state-variables $x_{i}(t)$, $i=1,\dots,n$. Just like in the RFM, $x_{i}(t)$ describes the normalized density at site $i$ at time $t$, where $x_{i}(t)=0$ [$x_{i}(t)=1$] means that the site is completely empty [full].

To state the dynamical equations describing the EFRBM we introduce more notation. Let $x_{0}(t)\equiv 1$, $x_{n+1}(t)\equiv 0$, and denote

[TABLE]

Then the EFRBM is described by

[TABLE]

where

[TABLE]

We now explain these equations. The term $g_{i}(x)$ represents the flow from site $i$ to site $i+1$, so Eq. (3) means that the change in the density at site $i$ is the inflow from site $i-1$ minus the outflow to site $i+1$. To explain Eq. (4), consider for example the case $i=2$ (and assume that $n\geq 4$). Then (4) yields

[TABLE]

The term $x_{2}$ means that the flow from site $2$ to site $3$ increases with the density at site $2$. The term $(1-x_{3})$ represents soft exclusion: as the density at site $3$ increases, the transition from site $2$ to site $3$ gradually decreases. The term $(1+(q-1)x_{4})$ represents the fact that the flow into site $3$ also depends on the density at site $4$: if $q>1$ [$q<1$] then the transition increases [decreases] with $x_{4}$, that is, the “particles” at site $4$ “attract” [“repel”] the particles that move from site $2$ to site $3$. The term $(1+(r-1)x_{1})$ is similar but represents an attachment/detachment force between the “particles” in sites $1$ and $2$.

Note that for $r=q=1$, $g_{i}(x)=\lambda_{i}x_{i}(1-x_{i+1})$, and thus in this case the EFRBM reduces to the RFM (see (1)). On the other hand, if $q=r=0$ then $g_{i}(x)=\lambda_{i}x_{i}(1-x_{i+1})(1-x_{i+2})(1-x_{i-1})$. This represents a kind of an “extended objects” RFM, as the transition from site $i$ to site $i+1$ decreases with the density in sites $i-1$, $i+1$, and $i+2$.

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

[TABLE]

where $\eta_{2}(t):=\lambda_{2}(1+(q-1)x_{4}(t))(1+(r-1)x_{1}(t))$. Note that this time-varying transition rate depends on $\lambda_{2}$ (i.e., the fixed site to site transition rate), and also on $r$ and $q$ 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 $x_{n}$ to the environment by

[TABLE]

This is the output rate at time $t$.

Example 1

The EFRBM with $n=3$ sites is given by:

[TABLE]

If $q=r=0$ then this becomes

[TABLE]

On the other-hand, for $q=1$ and $r=0$ (1) becomes

[TABLE]

and this system admits a continuum of steady-states, as $\begin{bmatrix}1&1&s\end{bmatrix}^{\prime}$ is a steady-state for all $s$.  $\square$

Following [38] (see also [39]), we view creating and breaking a pair of particles as opposite chemical transitions, so by detailed balance arguments: $\frac{q}{r}=\exp\left(\frac{E}{K_{B}T}\right)$, where $E$ is the interaction energy. As in [38], we also assume that $E$ is equally split between the creation and breaking processes, so

[TABLE]

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

[TABLE]

In this case, the EFRBM contains $n+2$ parameters: $\lambda_{0},\dots,\lambda_{n}$, and $r$ (as $q=1/r$). Note that if (11) holds then (4) becomes

[TABLE]

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 $x(t,a)$ denote the solution of the EFRBM at time $t$ for the initial condition $x(0)=a\in C^{n}$.

Invariance and persistence

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

Proposition 1

Assume that $q,r>0$. For any $\tau>0$ there exists $d=d(\tau)\in(0,1/2)$ such that

[TABLE]

for all $a\in C^{n}$, all $i\in\{1,\dots,n\}$, and all $t\geq 0$.

This means that all the trajectories of the EFRBM enter and remain in the interior of $C^{n}$ after an arbitrarily short time. In particular, both $C^{n}$ and $\operatorname{{\mathrm{I}nt}}(C^{n})$ 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 $[0,1]$ then this remains true for all time $t\geq 0$, so the equations “make sense” in this respect. Furthermore, after an arbitrarily short time the densities are all in $(0,1)$, 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:

[TABLE]

whose trajectories evolve on an invariant set $\Omega\subset\mathbb{R}^{n}$ that is compact and convex. Let $x(t,t_{0},a)$ denote the solution of (13) at time $t$ for the initial condition $x(t_{0})=a$. The dynamical system (13) is said to be contractive after a small overshoot (SO) [43] on $\Omega$ w.r.t. a norm $|\cdot|:\mathbb{R}^{n}\to\mathbb{R}_{+}$ if for any $\varepsilon>0$ there exists $\ell=\ell(\varepsilon)>0$ such that

[TABLE]

for all $a,b\in\Omega$ and all $t\geq t_{0}\geq 0$. Intuitively speaking, this means that any two trajectories of the system approach each other at an exponential rate $\ell$, but with an arbitrarily small overshoot of $1+\varepsilon$.

Let $|\cdot|_{1}:\mathbb{R}^{n}\to\mathbb{R}_{+}$ denote the $L_{1}$ norm, i.e. for $z\in\mathbb{R}^{n}$, $|z|_{1}=|z_{1}|+\dots+|z_{n}|$.

Proposition 2

The EFRBM with $q,r>0$ is SO on $C^{n}$ w.r.t. the $L_{1}$ norm, that is, for any $\varepsilon>0$ there exists $\ell=\ell(\varepsilon)>0$ such that

[TABLE]

for all $a,b\in C^{n}$ and all $t\geq 0$.

From a biological point of view this means the following. The state of the system at any time $t$ is a vector describing the density at each site at time $t$. We measure the distance between any two density vectors using the $L_{1}$ 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 $1_{n}$ [$0_{n}$] denote the column vector of $n$ ones [zeros].

Example 2

Consider the EFRBM with dimension $n=3$, and parameters $\lambda_{0}=1$, $\lambda_{1}=2$, $\lambda_{2}=3$, $\lambda_{3}=4$, $r=5$, and $q=1/5$. Fig. 3 depicts $|x(t,a)-x(t,b)|_{1}$, with $a=0_{3}$ and $b=1_{3}$, as a function of time for $t\in[0,2]$. It may be seen that the $L_{1}$ distance between the two trajectories goes to zero at an exponential rate.  $\square$

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 $\dot{x}=f(x)$. Since the compact and convex set $C^{n}$ is an invariant set of the dynamics, it contains at least one steady-state. That is, there exists $e=e(\lambda_{0},\dots,\lambda_{n},q,r)$ such that $f(e)=0_{n}$. By Proposition 1, $e\in\operatorname{{\mathrm{I}nt}}(C^{n})$. Using (14) with $b:=e$ yields the following result.

Corollary 1

Assume that $q,r>0$. Then the EFRBM admits a unique steady-state $e\in\operatorname{{\mathrm{I}nt}}(C^{n})$ that is globally asymptotically stable, i.e.

[TABLE]

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 $\lambda_{i}$, and the parameters $r$ and $q$, 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 $r,q>0$ cannot be dropped. Indeed, Eq. (1), corresponding to a EFRBM with $n=3$, $q=1$ and $r=0$, admits a continuum of steady-states.

Example 3

Fig. 4 depicts the trajectories of (3) with $n=3$, $\lambda_{0}=0.5$, $\lambda_{1}=0.8$, $\lambda_{2}=0.7$, $\lambda_{3}=0.6$, $r=1/2$, and $q=2$, for several initial conditions. It may be seen that all trajectories converge to a unique steady-state $e=\begin{bmatrix}0.8555&0.7881&0.4268\end{bmatrix}^{\prime}$. (All the numerical values in the simulations described in this paper are to four digit accuracy.)  $\square$

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.

Analysis of the steady-state

At steady-state, (i.e. for $x=e$) the left-hand side of all the equations in (3) is zero (i.e. $\dot{x}_{i}=0$, $i=1,\dots,n$), so $g_{i-1}(e)=g_{i}(e)$ for all $i$. This implies that

[TABLE]

and also that the steady-state flow satisfies

[TABLE]

In particular, $R=\lambda_{n}e_{n}(1+(r-1)e_{n-1})$ and since every $e_{i}\in(0,1)$, the steady-state flow is positive (i.e. a left-to-right flow) for any $r>0$.

Also, for the case $rq=1$ it follows from $R=\lambda_{0}(1-e_{1})(1+(q-1)e_{2})$ that for $r\geq 1$, $R\leq\lambda_{0}$, whereas for $r<1$ it follows from $R=\lambda_{n}e_{n}(1+(r-1)e_{n-1})$ that $R\leq\lambda_{n}$, so

[TABLE]

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 $\lambda_{i}$s by a parameter $c>0$ then $e$ will not change, i.e. $e(c\lambda)=e(\lambda)$. Thus, by (16) $R(c\lambda)=cR(\lambda)$, for all $c>0$, that is, the steady-state flow [density] is homogeneous of degree one [zero] w.r.t. the $\lambda_{i}$s.

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

Fact 2

Consider the EFRBM with $n=2$ and $q=1/r$. Define

[TABLE]

Then $e=\begin{bmatrix}e_{1}&e_{2}\end{bmatrix}^{\prime}$ is given by

[TABLE]

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

Let $\mathbb{R}^{n}_{++}$ denote the set of $n$ dimensional vectors with all entries positive. Let $v:=\begin{bmatrix}\lambda_{0}&\dots&\lambda_{n}&r&q\end{bmatrix}^{\prime}$ denote the set of parameters in the EFRBM with dimension $n$. The results above imply that there exists a function $h:\mathbb{R}^{n+3}_{++}\to\operatorname{{\mathrm{I}nt}}(C^{n})$ such that $e=h(v)$ is the unique steady-state of the EFRBM with parameters $v$.

Proposition 3

The function $h:\mathbb{R}^{n+3}_{++}\to\operatorname{{\mathrm{I}nt}}(C^{n})$ is analytic.

This result allows in particular to consider the derivatives of the steady-state density $e=e(v)$ and the steady-state output rate $R=R(v)$ w.r.t. small changes in some of the parameters $v$, 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 $r$ (and $q=1/r$) on the steady-state of the EFRBM.

Example 4

Consider a EFRBM with $n=2$ and rates $\lambda_{0}=\lambda_{1}=\lambda_{2}=1$. Fig. 5 depicts the steady-state output rate $R$ as a function of $r$. It may be seen that $R$ monotonically increases with $r$. In particular, for $r=1$ (i.e., the RFM) $R=0.3820$, wheres for $r=20$, $R=0.4778$, that is, the steady-state flow is increased by about $25\%$. When considering the comparison with the RFM, one should bear in mind that the EFRBM corresponds to an RFM with time-varying rates $\eta_{i}(t)$ that may effectively be much higher than the fixed rates $\lambda_{i}$. We assume that the energy that is needed to generate these higher rates comes from the additional interaction forces between the particles.  $\square$

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

Example 5

Consider the EFRBM with dimension $n=6$, $\lambda_{0}=1.0$, $\lambda_{1}=1.2$, $\lambda_{2}=0.9$, $\lambda_{3}=4.0$, $\lambda_{4}=0.2$, $\lambda_{5}=1.0$, and $\lambda_{6}=1.1$. Consider first the case $r=q=1$ (i.e., the RFM). The steady-state density is:

[TABLE]

and the corresponding steady-state flow is $R=0.1557$. Note that since $\lambda_{3}$ is high and $\lambda_{4}$ is low, $e_{1},e_{2},e_{3},e_{4}\gg e_{5},e_{6}$, indicating a traffic jam at site $4$. Consider now the case $r=5$ (i.e. $q=1/5$). The steady-state density is now

[TABLE]

and the corresponding steady-state output is $\tilde{R}=0.2820$. Note that now the density at site $4$ decreased relative to the $r=1$ case, and that $\tilde{R}>R$. Note also that $\sum_{i=1}^{6}e_{i}=3.7602>\sum_{i=1}^{6}\tilde{e}_{i}=2.4403$. This means that the introduction of a “neighbor-repelling” force (i.e. $r>1$) 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 $r\in[1,10]$. It may be observed that $e_{i}$, $i=1,\dots,5$, monotonically decreases with $r$, and that $e_{6}$ slightly increases with $r$. Note that since the occupancy at site $6$ is not affected by $q$, but only by $r$, increasing $r$ should indeed increase $e_{6}$.  $\square$

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 $r\to 0$ (so $q=\frac{1}{r}\to\infty$) and $r\to\infty$ (so $q=\frac{1}{r}\to 0$).

The case $r\to 0$

Intuitively speaking, a low value of $r$ corresponds to: (1) a strong attachment between existing nearest neighbors (small $r$); and (2) a high tendency for moving forward if this involves creating new neighbors (large $q$). 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 $n=6$ and rates $\lambda_{i}=1$, $i=0,\dots,6$. For $r=0.1$ (recall that $q=1/r$), the steady-state values are:

[TABLE]

For $r=0.01$, the steady-state values are:

[TABLE]

For $r=0.005$, the steady-state density values are:

[TABLE]

Fig. 7 depicts the steady-state values for the three $r$ values. It may be observed that as $r$ 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.  $\square$

We now rigorously analyze the case $r\to 0$ for the EFRBM with $n=2$ and $n=3$.

Example 7

Consider the EFRBM with $n=2$ and $q=1/r$. Expanding $e_{2}$ and $e_{1}$ in (2) as a Taylor series in $r$ yields

[TABLE]

where every $o(r)$ denotes a function $f(r)$ satisfying $\lim_{r\to 0}\frac{f(r)}{r}=0$. Thus, $R=\lambda_{1}e_{1}(1-e_{2})=\lambda_{2}r+o(r)$. This implies in particular that

[TABLE]

Thus, when $r\to 0$, both steady-state densities go to one,111Note that (19) implies that $e_{1}$ goes to one faster than $e_{2}$. that is, the sites become completely full, and consequently the steady-state output rate goes to zero.  $\square$

The next result analyzes the case $n=3$.

Proposition 4

The steady-state densities in the EFRBM with $n=3$ satisfy

[TABLE]

and

[TABLE]

Note that this implies that

[TABLE]

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

Summarizing, as $r$ goes to [math] 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, $r\to\infty$.

The case $r\to\infty$

A large value of $r$ corresponds to: (1) strong repulsion between existing nearest neighbors (large $r$); and (2) a low tendency for moving forward if this involves creating new neighbors (small $q$). 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 $n=6$ sites and rates $\lambda=\begin{bmatrix}1&1.2&0.8&0.95&1.1&0.75&1.15\end{bmatrix}^{\prime}.$ For $r=1$ (recall that $q=1/r$),

[TABLE]

For $r=1,000$,

[TABLE]

For $r=10,000$,

[TABLE]

Fig. 8 depicts these steady-state values for the three $r$ values. Note that the steady-state values for $r=1,000$ and $r=10,000$ cannot be distinguished. It may be observed that the values $e_{j}(r)$, $j=2,4,6$, decrease to zero as $r$ 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 $R$ as $r$ increases.  $\square$

We now rigorously analyze the case $r\to\infty$ for the EFRBM with $n=2$ and $n=3$.

Example 9

Consider the EFRBM with $n=2$. Expanding $e$ in (2) as a Taylor series in $q=1/r$ yields

[TABLE]

so

[TABLE]

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

Proposition 5

The steady-state densities in the EFRBM with $n=3$ satisfy

[TABLE]

with $a_{1},a_{3}\in(0,1)$, and

[TABLE]

Note that this implies that

[TABLE]

so again as $r\to\infty$ the density at site $2$ 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 $\lambda_{i}$ between sites. The nearest-neighbor interactions between the particles are modeled using two parameters: $r$ that represents the tendency of a moving particle to break from an existing neighbor, and $q$ 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 $t$ represent the probability of each configuration at time $t$, and the number of possible configurations grows exponentially with the number of sites $n$. On the other-hand, the EFRBM includes $n$ (nonlinear) ODEs for $n$ 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 $r$ is very large and $q$ 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 $r$ (or $q$) means very high effective transition rates, but our results suggest that even for $r$ 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 $C^{n}$ is an invariant set of the dynamics follows immediately from the equations of the EFRBM. Let

[TABLE]

with the $z_{i}$s defined in (2). By (3), the EFRBM can be written as

[TABLE]

This is just the RFM (see (1)), but with time-varying rates $\eta_{i}(t)$. Let $a_{i}:=\min\{1,q\}\min\{1,r\}\lambda_{i}$, and $b_{i}:=\max\{1,q\}\max\{1,r\}\lambda_{i}$. It follows from (24) that $a_{i}\leq\eta_{i}(t)\leq b_{i}$ for all $i$ and for all $t\geq 0$. Note that for $r,q>0$ every $a_{i}$ 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]. $\square$

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 $C^{n}$. Also, Proposition 4 in [43] implies that for the EFRBM the properties of SOST and SO are equivalent, and this completes the proof. $\square$

Proof of Fact 2. Consider the EFRBM with $n=2$ and $q=1/r$. Then (Analysis of the steady-state) becomes

[TABLE]

This yields

[TABLE]

and

[TABLE]

with $a_{1}$ defined in (17) and $a_{2}:=-\lambda_{1}-\lambda_{2}-a_{1}$. The feasible solution (i.e. the one satisfying $e_{1},e_{2}\in(0,1)$ for any set of parameter values) is given by

[TABLE]

and (26). $\square$

Proof of Prop. 3. To emphasize the dependence on the parameters, write the EFRBM as $\dot{x}=f(x;v)$, where $v:=\begin{bmatrix}\lambda_{0}&\dots&\lambda_{n}&r&q\end{bmatrix}^{\prime}$. Note that $f$ is an analytic function. Then the steady-state satisfies the relation $f(e;v)=0$. The Jacobian matrix of this relation with respect to $x$ is

[TABLE]

which is just the Jacobian of the dynamics. Fix $v^{0}\in\mathbb{R}^{n+3}_{++}$ and let $e^{0}\in\operatorname{{\mathrm{I}nt}}(C^{n})$ denote the corresponding steady-state, that is, $f(e^{0};v^{0})=0$ and $e^{0}=h(v^{0})$. Suppose that there exists a matrix measure $\mu:\mathbb{R}^{n\times n}\to\mathbb{R}$ such that $\mu(J(e^{0};v^{0}))<0$. This implies in particular that $J(e^{0},v^{0})$ is Hurwitz (see e.g. [46]), so it is not singular and invoking the implicit function theorem implies that the mapping $h$ is analytic. It follows from the results in [28] that such a matrix measure $\mu$ indeed exists, and this completes the proof. $\square$

Proof of Prop. 4. Expand $e_{i}$, $i=1,2,3$, as

[TABLE]

Recall that the steady-state equations are given by $R(e)=g_{0}(e)=g_{1}(e)=\dots=g_{3}(e)$, with the $g_{i}$s given in (4). Substituting (27) yields

[TABLE]

Since $R(e)$ is bounded, we conclude that

[TABLE]

Assume for the moment that $a_{2}=0$. Then

[TABLE]

and this implies that $a_{3}=0$. Now, $g_{1}(e)=\lambda_{1}(1+b_{3})+\dots$ and combining this with (29) yields $b_{3}=-1$. Thus, e_{3}=a_{3}+b_{3}r+c_{3}r^{2}+\operatorname{o}(r^{2})=-r+\operatorname{o}\bigl{(}r\bigr{)}, and this is a contradiction as $e_{3}(r)$ will be strictly negative for any $r>0$ sufficiently small. We conclude that $a_{2}\not=0$, so (28) yields $a_{1}=1$, and also $(1-a_{2})a_{3}=0$. Suppose that $a_{3}=0$. Then

[TABLE]

It follows that $a_{2}b_{1}=(1+b_{3})(1-a_{2})=0$. Since we already know that $a_{2}\not=0$, $b_{1}=0$. The case $b_{3}=-1$ is impossible, as then $e_{3}(r)<0$ for $r>0$ sufficiently small, so $a_{2}=1$. But then $R(e)=g_{2}(e)=-\lambda_{2}r+\dots$ and this is a contradiction. We conclude that $a_{3}\not=0$, so (28) yields $a_{2}=1$. Summarizing, we have $a_{1}=a_{2}=1$. Now,

[TABLE]

This gives $b_{1}=0$ and $b_{2}a_{3}=0$. Since we already know that $a_{3}\not=0$, $b_{2}=0$. Now,

[TABLE]

Equating the coefficients here yields $a_{3}=\frac{\lambda_{2}}{\lambda_{2}+\lambda_{3}}$, $c_{1}=\frac{-\lambda_{2}\lambda_{3}}{\lambda_{0}(\lambda_{2}+\lambda_{3})}$, and $c_{2}=-\lambda_{3}/\lambda_{1}$. Since we know that the steady-state equations admit a unique solution this yields (4), and the equation $R(e)=g_{0}(e)$ yields (21). $\square$

Proof of Prop. 5. Expand $e_{i}$, $i=1,2,3$, as

[TABLE]

Recall that the steady-state equations are given by $R(e)=g_{0}(e)=g_{1}(e)=\dots=g_{3}(e)$, with the $g_{i}$s given in (4). Substituting (30) yields

[TABLE]

This implies that

[TABLE]

Assume for the moment that $a_{2}\not=0$. Then $a_{3}=a_{1}=0$. This yields

[TABLE]

This implies that $a_{2}=1$ and $b_{1}=-1$. This yields e_{1}(r)=-r+\operatorname{o}\bigl{(}r\bigr{)} which is a contradiction.

We conclude that $a_{2}=0$. Now,

[TABLE]

Equating the coefficients here yields the following. First, $a_{1}\not=0$, and since $e_{1}(q)=a_{1}+\dots$, this implies that $a_{1}>0$. Second, if $a_{3}=1$ then $a_{1}=1$ and $b_{2}=-1$ which is a contradiction as then $e_{2}(q)=-q+o(q)$. Thus, $a_{3}\not=1$ and this implies that $a_{1}\not=1$. We conclude that $a_{1},a_{3}\in(0,1)$. Now the equations for $g_{1}$ and $g_{2}$ yield $b_{2}=\lambda_{1}/\lambda_{2}$. This proves (5). Expanding $g_{0}$ up to order one in $q$, and using $R=g_{0}$ yields (23). $\square$