A Deterministic Model for One-Dimensional Excluded Flow with Local Interactions
AA Deterministic Model for One-DimensionalExcluded Flow with Local Interactions
Yoram Zarai , Michael Margaliot , ∗ , Anatoly B. Kolomeisky ∗ Corresponding Author. Email: [email protected]
Abstract
Natural phenomena frequently involve a very large number of interacting molecules mov-ing in confined regions of space. Cellular transport by motor proteins is an example ofsuch collective behavior. We derive a deterministic compartmental model for the unidi-rectional flow of particles along a one-dimensional lattice of sites with nearest-neighborinteractions between the particles. The flow between consecutive sites is governed bya “soft” simple exclusion principle and by attracting or repelling forces between neigh-boring particles. Using tools from contraction theory, we prove that the model admitsa unique steady-state and that every trajectory converges to this steady-state. Analysisand simulations of the effect of the attracting and repelling forces on this steady-statehighlight 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 dynamicprocesses.
Introduction
Biological processes are governed by complex interactions between multiple particlesthat are confined in special compartments [1]. One of the most important examples ofsuch 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 asbiological molecular motors, can catalyze the reaction of adenosine triphosphate (ATP)hydrolysis, while at the same time converting the energy produced during this chemicalreaction 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 largegroups, suggesting that the interactions between the motors cannot be ignored [3,4]. Un-derstanding the collective behavior of molecular motors is critical for uncovering mech-anisms of complex biological processes [2, 5, 6]. From a theoretical point of view, intra-cellular transport processes are usually described using non-equilibrium multi-particlelattice models [3]. In these models, the molecular motors are typically represented byparticles that hop along the lattice, and the lattice sites model the binding locations ofthe motors along the filaments (or tracks). For a general review on transport and trafficphenomena in biological systems see for example [2–5]. a r X i v : . [ q - b i o . S C ] M a y standard model from non-equilibrium statistical mechanics for molecular motorstraffic (and numerous other processes) is the totally asymmetric simple exclusion pro-cess (TASEP) [7–9]. This is also the standard model for ribosome flow during mRNAtranslation (see, e.g. [8, 10, 11]). In TASEP, particles hop randomly along a unidirec-tionally ordered lattice of sites. Simple exclusion means that a particle cannot moveinto a site that is already occupied by another particle, and thus each site can be eitherempty or occupied by a single particle. This models moving biological particles likeribosomes and motor proteins that have volume and thus cannot overtake a movingparticle in front of them. This hard exclusion principle creates an intricate indirectcoupling between the particles. In particular, a slowly moving particle may lead to theformation of a traffic jam behind it.To describe moving biological molecules with large sizes, a version of TASEP withextended objects has been introduced and analyzed [12–14]. In this model, each particlecovers (cid:96) > i, . . . , i + (cid:96) − i , andit can hop to site i + 1 provided that site i + (cid:96) is empty. This is used, for example, formodeling mRNA translation as it is known that every ribosome (the particle) coversseveral codons (sites) along the mRNA molecule [13].There exist two versions of TASEP that differ by their boundary conditions. In TASEPwith open boundary conditions the two sides of the chain are connected to two particlereservoirs 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 openboundary homogeneous TASEP (HTASEP), all the transition rates within the latticeare assumed to be equal and normalized to one, and thus the model is specified by aninput rate α , an exit rate β , and a parameter N denoting the number of sites along thelattice. In TASEP with periodic boundary conditions the chain is closed into a ring, anda particle that hops from the last site returns to the first site. TASEP has been widelyutilized for studying various natural and artificial processes, including vehicular trafficflow, mRNA translation, surface growth, communication networks, and more [3, 15].Ref. [16] used HTASEP with periodic boundary conditions to analyze transport ona lattice in the presence of local interactions between particles and substrate, illustrat-ing the effect of local conformation of the substrate on the characteristics of the flowof molecular motors. TASEP with particle interactions and with periodic boundaryconditions was studied in [17], and with open boundary conditions in [18–21]. Specif-ically, the authors in [20, 21] proposed a modified TASEP model that incorporates therealistic observed feature of nearest-neighbor interactions. In this model, the transitionrate in every site along the lattice depends on the states of four consecutive sites. Theirconclusions were that weak repulsive interaction results in maximal flux, and that themolecular motors are influenced more strongly by attractive interactions.Unfortunately, rigorous analysis of TASEP is non-trivial, and exact solutions existonly in special cases, for example when considering the model with the homogeneous rates (HTASEP). Typically, the non-homogeneous case and cases that include other localinteractions are only studied via various approximations and extensive Monte Carlocomputer simulations. These simulations are run until convergence to a (stochastic)steady-state, yet without a rigorous proof that convergence indeed takes place for allthe feasible parameter values.In this paper, we introduce a new deterministic model for the flow of motor proteinsalong a one-dimensional lattice of sites with nearest-neighbor interactions between themotors. The flow of the motor proteins is unidirectional, and it satisfies a “soft” simpleexclusion principle. The nearest-neighbor effect is modeled by two “force” interactionswith parameters q and r . It is more convenient to explain the effect of these interactionsin “particle-like” terms, although in the new model the density in every site takes valuesin the range [0 ,
1] (and not { , } ).onsider a transition of a particle from site i to site i + 1. If site i + 2 is alreadyoccupied then the rate of movement depends on a parameter q ≥ new neighbors. A value q > q < i + 2. On the other-hand, if site i − r ≥ old neighbors. A value r > r <
1] means that the particle will tend [not] to hop forward, as there is a strongrepulsion [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 [whenbreaking from old neighbors].An important advantage of our model is that it is highly amenable to rigorousanalysis even for non-homogenous transition rates. We prove, for example, that thedynamics always converges to a steady-state density along the lattice. Thus, the flowalso converges to a steady-state value. This steady-state depends on the lattice size, thetransition 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 attractingand repelling forces on this steady-state highlight the crucial role that these forces mayplay in increasing the steady-state flow, and reveal that this increase stems from thealleviation of traffic jams along the lattice. It is well-known that molecular motorsindeed form traffic jams and that these have important biological implications (see,e.g. [22–24]). In particular, analysis and simulations of the model reveal a new regimethat may be interpreted as the “opposite” of a traffic jam along the lattice.Our approach extends a deterministic mathematical model that has been used fordescribing and analyzing the flow of ribosomes along the mRNA molecule during theprocess 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 unidirec-tional 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 condi-tions [3, Section 4.9.7] [7, p. R345]. The RFM includes n +1 parameters: λ > λ n > λ i > i = 1 , . . . , n −
1, the transition ratefrom site i to site i + 1. The state variable x i ( t ) : R + → [0 , i = 1 , . . . , n , describesthe 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 ) := (cid:2) x ( t ) . . . x n ( t ) (cid:3) (cid:48) describes the density profile along the chain attime t . The output rate at time t is R ( t ) := λ n x n ( t ) (see Fig. 1).Let x ( t ) ≡
1, and x n +1 ( t ) ≡
0. The dynamics of the RFM with n sites is given bythe following set of n nonlinear ODEs:˙ x i = λ i − x i − (1 − x i ) − λ i x i (1 − x i +1 ) , i = 1 , . . . , n. (1)This can be explained as follows. The flow of material from site i to site i + 1 at time t is λ i x i ( t )(1 − x i +1 ( t )). This flow increases with the density at site i , and decreases assite i + 1 becomes fuller. This corresponds to a “soft” version of a simple exclusionprinciple. Note that the maximal possible flow from site i to site i + 1 is the transitionrate λ i . Thus Eq. (1) simply states that the change in the density at site i at time t isthe 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 C n := { x ∈ R n : x i ∈ [0 , , i = 1 , . . . , n } . igure 1. The RFM models unidirectional flow along a chain of n sites. The statevariable x i ( t ) ∈ [0 ,
1] represents the density at site i at time t . The parameter λ i > i to site i + 1, with λ > λ n >
0] controllingthe initiation [exit] rate. The output rate at time t is R ( t ) := λ n x n ( t ).Let Int( C n ) [ ∂C n ] denote the interior [boundary] of C n . Ref. [26] has shown that theRFM is a tridiagonal cooperative dynamical system [27], and consequently (1) admits a unique steady-state density e = e ( λ , . . . , λ n ) ∈ Int( C n ) that is globally asymptoticallystable, that is, lim t →∞ x ( t, a ) = e for all a ∈ C n (see also [28]). This means thattrajectories 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 := λ n e n .An important advantage of the RFM (e.g. as compared to TASEP) is that it isamenable to mathematical analysis using tools from systems and control theory. Fur-thermore, most of the analysis hold for the general, non-homogeneous case (i.e. the casewhere the transition rates λ i differ from one another). For more on the analysis of theRFM and its biological implications, see [26, 28–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 thestate of its neighboring sites. A parameter r [ q ] controls the binding/repelling forcesbetween two existing [new] neighbors. We refer to the new model as the excluded flowwith local repelling and binding model (EFRBM). It is important to note that this issignificantly 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 byits two nearest neighbors sites, whereas in the EFRBM the dynamics is directly affectedby the density in four neighboring sites. Thus, unlike the RFM, the EFRBM is not atridiagonal system. Also, the RFM has been used to model ribosome flow, whereas herewe apply the EFRBM to study the flow of motor proteins.We show that the EFRBM is a contractive dynamical system. This holds for anyset 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 correspondsto a unique steady-state output rate. We analyze the behavior of this steady-stateunder 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 alarge q = 1 /r ) leads to a small output rate. Analysis and simulations show that this isdue 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 goesto zero. This “separation of densities” is the “opposite” of a traffic jam. These resultshighlight the impact of traffic jams on the output rate.The remainder of this paper is organized as follows. The next section describesthe EFRBM. The following two sections describe our main analysis results and their bi-ological implications. This includes analysis of the asymptotic behavior of the EFRBM, igure 2. A schematic explanation of the transition flow from site i to site i + 1 inthe EFRBM. Upper-left: when both sites i − i + 2 do not contain particles, thetransition rate is λ i . Lower-left: when site i − i + 2 does, the transition rate is λ i q . Upper-right: when site i − i + 2 does not, the transition rate is λ i r . Lower-right: when bothsites i − i + 2 contain particles, the transition rate is λ i rq .and the effects of the nearest-neighbor interactions on the steady-state behavior ofthe EFRBM. The final section summarizes and describes several directions for furtherresearch. To increase the readability of this paper, all the proofs are placed in theAppendix. The EFRBM
The EFRBM with n sites includes n + 3 parameters: • λ i > i = 0 , . . . , n , controls the transition rate from site i to site i + 1, where λ [ λ n ] controls the input [output] rate. • r ≥ • q ≥ i to site i + 1, andthe 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 inthe range [0 ,
1] and not a binary choice { , } like in TASEP. If both sites i − i + 2do not contain particles, the transition rate is simply λ i , as in the RFM. If a particle islocated at site i − i + 2] but site i + 2 [ i −
1] is empty then the transition rate is λ i r [ λ i q ]. If both sites contain particles the transition rate is λ i rq .The EFRBM also includes n state-variables x i ( t ), i = 1 , . . . , 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 nota-tion. Let x ( t ) ≡ x n +1 ( t ) ≡
0, and denote z i ( t ) := (cid:40) x i ( t ) , i = 1 , . . . , n, , otherwise . (2)hen the EFRBM is described by˙ x i = g i − ( x ) − g i ( x ) , i = 1 , . . . , n, (3)where g i ( x ) := λ i x i (1 − x i +1 )(1 + ( q − z i +2 )(1 + ( r − z i − ) . (4)We now explain these equations. The term g i ( x ) represents the flow from site i tosite i + 1, so Eq. (3) means that the change in the density at site i is the inflow fromsite i − i + 1. To explain Eq. (4), consider for example thecase i = 2 (and assume that n ≥ g ( x ) = λ x (1 − x )(1 + ( q − x )(1 + ( r − x ) . (5)The term x means that the flow from site 2 to site 3 increases with the density atsite 2. The term (1 − x ) represents soft exclusion: as the density at site 3 increases, thetransition from site 2 to site 3 gradually decreases. The term (1 + ( q − x ) representsthe fact that the flow into site 3 also depends on the density at site 4: if q > q <
1] thenthe transition increases [decreases] with x , that is, the “particles” at site 4 “attract”[“repel”] the particles that move from site 2 to site 3. The term (1 + ( r − x ) is similarbut represents an attachment/detachment force between the “particles” in sites 1 and 2.Note that for r = q = 1, g i ( x ) = λ i x i (1 − x i +1 ), and thus in this case the EFRBMreduces to the RFM (see (1)). On the other hand, if q = r = 0 then g i ( x ) = λ i x i (1 − x i +1 )(1 − x i +2 )(1 − x i − ). This represents a kind of an “extended objects” RFM, asthe transition from site i to site i + 1 decreases with the density in sites i − i + 1,and i + 2. Remark 1
It is useful to think of the EFRBM as an RFM with time-varying transitionrates. For example, we can write (5) as g ( x ( t )) = η ( t ) x ( t )(1 − x ( t )) , where η ( t ) := λ (1 + ( q − x ( t ))(1 + ( r − x ( t )). Note that this time-varyingtransition rate depends on λ (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 theinteraction forces between the moving particles.We denote the flow from site x n to the environment by R ( t ) := λ n x n ( t )(1 + ( r − x n − ( t )) . (6)This is the output rate at time t . Example 1
The EFRBM with n = 3 sites is given by:˙ x = λ (1 − x )(1 + ( q − x ) − λ x (1 − x )(1 + ( q − x ) , ˙ x = λ x (1 − x )(1 + ( q − x ) − λ x (1 − x )(1 + ( r − x ) , ˙ x = λ x (1 − x )(1 + ( r − x ) − λ x (1 + ( r − x ) . (7)If q = r = 0 then this becomes˙ x = λ (1 − x )(1 − x ) − λ x (1 − x )(1 − x ) , ˙ x = λ x (1 − x )(1 − x ) − λ (1 − x ) x (1 − x ) , ˙ x = λ (1 − x ) x (1 − x ) − λ (1 − x ) x . (8)n the other-hand, for q = 1 and r = 0 (7) becomes˙ x = λ (1 − x ) − λ x (1 − x ) , ˙ x = λ x (1 − x ) − λ (1 − x ) x (1 − x ) , ˙ x = λ (1 − x ) x (1 − x ) − λ (1 − x ) x , (9)and this system admits a continuum of steady-states, as (cid:2) s (cid:3) (cid:48) is a steady-state forall s . (cid:3) Following [38] (see also [39]), we view creating and breaking a pair of particles as op-posite chemical transitions, so by detailed balance arguments: qr = exp (cid:16) EK B T (cid:17) , where E is the interaction energy. As in [38], we also assume that E is equally split between thecreation and breaking processes, so q = exp (cid:18) E K B T (cid:19) , r = exp (cid:18) − E K B T (cid:19) . (10)This has a clear physical meaning. If E > q >
1) since the energy of the system decreasesby E . On the other-hand, breaking out of the cluster increases the energy by E and thetransition rate is thus slowed down ( r < E < q < r >
1. Note that (10) implies in particularthat rq = 1 . (11)In this case, the EFRBM contains n + 2 parameters: λ , . . . , λ n , and r (as q = 1 /r ).Note that if (11) holds then (4) becomes g i ( x ) = λ i x i (1 − x i +1 )(1 − r − r z i +2 )(1 + ( r − z i − ) . (12)The next section derives several theoretical results on the dynamical behavior ofthe 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 ∈ C n . Invariance and persistence
The next result shows that the n -dimensional unit cube C n is an invariant set ofthe 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 boundaryof C n “immediately enters” C n . This is a technical result, but it is important as in theinterior of C n the EFRBM admits several useful properties. Proposition 1
Assume that q, r > . For any τ > there exists d = d ( τ ) ∈ (0 , / such that d ≤ x i ( t + τ, a ) ≤ − d, for all a ∈ C n , all i ∈ { , . . . , n } , and all t ≥ . his means that all the trajectories of the EFRBM enter and remain in the interiorof C n after an arbitrarily short time. In particular, both C n and Int( C n ) are invariantsets of the EFRBM dynamics.From a biological point of view this means that if the system is initiated such thatevery density is in [0 ,
1] then this remains true for all time t ≥
0, so the equations“make sense” in this respect. Furthermore, after an arbitrarily short time the densitiesare all in (0 , Contraction
Differential analysis and in particular contraction theory proved to be a powerful toolfor analyzing the asymptotic behavior of nonlinear dynamical systems. In a contractivesystem, trajectories that emanate from different initial conditions approach each otherat an exponential rate [40–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 dynamicalsystem: ˙ x ( t ) = f ( t, x ( t )) , (13)whose trajectories evolve on an invariant set Ω ⊂ R n that is compact and convex.Let x ( t, t , a ) denote the solution of (13) at time t for the initial condition x ( t ) = a .The dynamical system (13) is said to be contractive after a small overshoot (SO) [43]on Ω w.r.t. a norm | · | : R n → R + if for any ε > (cid:96) = (cid:96) ( ε ) > | x ( t, t , a ) − x ( t, t , b ) | ≤ (1 + ε ) exp( − ( t − t ) (cid:96) ) | a − b | , for all a, b ∈ Ω and all t ≥ t ≥
0. Intuitively speaking, this means that any twotrajectories of the system approach each other at an exponential rate (cid:96) , but with anarbitrarily small overshoot of 1 + ε .Let | · | : R n → R + denote the L norm, i.e. for z ∈ R n , | z | = | z | + · · · + | z n | . Proposition 2
The EFRBM with q, r > is SO on C n w.r.t. the L norm, that is,for any ε > there exists (cid:96) = (cid:96) ( ε ) > such that | x ( t, a ) − x ( t, b ) | ≤ (1 + ε ) exp( − (cid:96)t ) | a − b | , (14) for all a, b ∈ C n and all t ≥ . From a biological point of view this means the following. The state of the systemat any time t is a vector describing the density at each site at time t . We measure thedistance between any two density vectors using the L vector norm. Suppose that weinitiate the system with two different densities. This generates two different solutionsof the dynamical system. The distance between these solutions decreases with time atan exponential rate.The next example demonstrates this contraction property. Let 1 n [0 n ] denote thecolumn vector of n ones [zeros]. Example 2
Consider the EFRBM with dimension n = 3, and parameters λ = 1, λ = 2, λ = 3, λ = 4, r = 5, and q = 1 /
5. Fig. 3 depicts | x ( t, a ) − x ( t, b ) | , with a = 0 and b = 1 , as a function of time for t ∈ [0 , L distance between the two trajectories goes to zero at an exponential rate. (cid:3) Prop. 2 implies that the EFRBM satisfies several important asymptotic properties.These are described in the following subsections. igure 3.
The distance | x ( t, a ) − x ( t, b ) | as a function of time for the EFRBM inExample 2. Global asymptotic stability
Write the EFRBM (3) as ˙ x = f ( x ). Since the compact and convex set C n is aninvariant set of the dynamics, it contains at least one steady-state. That is, thereexists e = e ( λ , . . . , λ n , q, r ) such that f ( e ) = 0 n . By Proposition 1, e ∈ Int( C n ).Using (14) with b := e yields the following result. Corollary 1
Assume that q, r > . Then the EFRBM admits a unique steady-state e ∈ Int( C n ) that is globally asymptotically stable, i.e. lim t →∞ x ( t, a ) = e, for all a ∈ C n . This means that any solution of the EFRBM converges to a unique steady-statedensity (and thus a unique steady-state output rate) that depends on the rates λ i ,and the parameters r and q , but not on the initial condition. From a biological pointof view, this means that the system always converges to a steady-state density and acorresponding steady-state output rate, and thus it makes sense to study how thesedepend on the various parameters.Note that the assumption that r, q > 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, λ = 0 . λ = 0 . λ = 0 . λ = 0 . r = 1 /
2, and q = 2, for several initial conditions. It may be seenthat all trajectories converge to a unique steady-state e = (cid:2) . . . (cid:3) (cid:48) .(All the numerical values in the simulations described in this paper are to four digitaccuracy.) (cid:3) The rigorous proof that every trajectory converges to a steady-state is important, asit 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-stateoutput rate, and explore how these are related to the various parameters of the model.
Figure 4.
Trajectories of the EFRBM in Example 3 for seven arbitrary initialconditions. The steady-state e is denoted by ∗ . 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.˙ x i = 0, i = 1 , . . . , n ), so g i − ( e ) = g i ( e ) for all i . This implies that λ (1 − e )(1 + ( q − e )= λ e (1 − e )(1 + ( q − e )= λ e (1 − e )(1 + ( q − e )(1 + ( r − e )= λ e (1 − e )(1 + ( q − e )(1 + ( r − e )... (15)= λ n − e n − (1 − e n )(1 + ( r − e n − )= λ n e n (1 + ( r − e n − ) , and also that the steady-state flow satisfies R = g i ( e ) , i = 0 , . . . , n. (16)In particular, R = λ n e n (1 + ( r − e n − ) and since every e i ∈ (0 , r > rq = 1 it follows from R = λ (1 − e )(1 + ( q − e ) that for r ≥ R ≤ λ , whereas for r < R = λ n e n (1 + ( r − e n − ) that R ≤ λ n , so R ≤ max { λ , λ n } . This means in particular that the output rate is always bounded.
Fact 1
It follows from (15) that if we multiply all the λ i s by a parameter c > then e will not change, i.e. e ( cλ ) = e ( λ ) . Thus, by (16) R ( cλ ) = cR ( λ ) , for all c > , that is,the steady-state flow [density] is homogeneous of degree one [zero] w.r.t. the λ i s. In the spacial case n = 2 the steady-state equations (15) can be solved in closed-form. act 2 Consider the EFRBM with n = 2 and q = 1 /r . Define a := (1 − r )( λ r + λ ) + λ λ λ . (17) Then e = (cid:2) e e (cid:3) (cid:48) is given by e = λ + λ + a − (cid:112) ( λ + λ + a ) − a λ a ,e = λ e λ + ( λ (1 − r ) − λ ) e . (18) Note that even in this case the expression for e is non-trivial. Let R n ++ denote the set of n dimensional vectors with all entries positive. Let v := (cid:2) λ . . . λ n r q (cid:3) (cid:48) denote the set of parameters in the EFRBM with dimension n .The results above imply that there exists a function h : R n +3++ → Int( C n ) such that e = h ( v ) is the unique steady-state of the EFRBM with parameters v . Proposition 3
The function h : R n +3++ → Int( 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 theparameters v , that is, the sensitivity of the steady-state w.r.t. small changes in theparameters. 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 λ = λ = λ = 1. Fig. 5 depictsthe steady-state output rate R as a function of r . It may be seen that R monotonicallyincreases with r . In particular, for r = 1 (i.e., the RFM) R = 0 . r = 20, R = 0 . η i ( t ) that may effectively be much higher than thefixed rates λ i . We assume that the energy that is needed to generate these higher ratescomes from the additional interaction forces between the particles. (cid:3) The next example demonstrates that the increase in R as r increases is because theneighbor-repelling forces lead to an alleviation of traffic jams. Example 5
Consider the EFRBM with dimension n = 6, λ = 1 . λ = 1 . λ = 0 . λ = 4 . λ = 0 . λ = 1 .
0, and λ = 1 .
1. Consider first the case r = q = 1 (i.e., theRFM). The steady-state density is: e = (cid:2) . . . . . . (cid:3) (cid:48) , and the corresponding steady-state flow is R = 0 . λ is high and λ is low, e , e , e , e (cid:29) e , e , indicating a traffic jam at site 4. Consider now the case r = 5 (i.e. q = 1 / e = (cid:2) . . . . . . (cid:3) (cid:48) , igure 5. Steady-state output rate R as a function of r ∈ [0 . ,
20] for a EFRBMwith n = 2, λ i = 1 for all i , and q = 1 /r . Note that the value for r = 1 is thesteady-state output rate in the RFM.and the corresponding steady-state output is ˜ R = 0 . r = 1 case, and that ˜ R > R . Note also that (cid:80) i =1 e i =3 . > (cid:80) i =1 ˜ e i = 2 . 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 ∈ [1 , e i , i = 1 , . . . ,
5, monotonically decreases with r , and that e 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 . (cid:3) 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 → q = r → ∞ ) and r → ∞ (so q = r → The case r → r corresponds to: (1) a strong attachment betweenexisting nearest neighbors (small r ); and (2) a high tendency for moving forward if thisinvolves creating new neighbors (large q ). As we will see this leads to the formation oftraffic jams and, consequently, to a sharp decrease in the output rate. Example 6
Consider a EFRBM with dimension n = 6 and rates λ i = 1, i = 0 , . . . , r = 0 . q = 1 /r ), the steady-state values are: e = (cid:2) . . . . . . (cid:3) (cid:48) , R = 0 . . For r = 0 .
01, the steady-state values are: e = (cid:2) . . . . . . (cid:3) (cid:48) , R = 0 . . r e e e e e e Figure 6.
Steady-state densities e i as a function of r ∈ [1 ,
10] for a EFRBMwith n = 6, λ = 1 . λ = 1 . λ = 0 . λ = 4 . λ = 0 . λ = 1 . λ = 1 . q = 1 /r . Note that as r increases all densities become much smaller than one,that is, there are no traffic jams.For r = 0 . e = (cid:2) . . . . . . (cid:3) (cid:48) , R = 0 . . Fig. 7 depicts the steady-state values for the three r values. It may be observed thatas r decreases the density in the first five sites increases to one, i.e. these sites becomecompletely full, and the output rate goes to zero. Note that this highlights the negativeeffect of traffic jams on the output rate. (cid:3) We now rigorously analyze the case r → n = 2 and n = 3. Example 7
Consider the EFRBM with n = 2 and q = 1 /r . Expanding e and e in (18) as a Taylor series in r yields e = 1 − λ λ r + o ( r ) , e = 1 + o ( r ) , (19)where every o ( r ) denotes a function f ( r ) satisfying lim r → f ( r ) r = 0. Thus, R = λ e (1 − e ) = λ r + o ( r ). This implies in particular thatlim r → e = lim r → e = 1 , lim r → R = 0 . Thus, when r →
0, both steady-state densities go to one, that is, the sites becomecompletely full, and consequently the steady-state output rate goes to zero. (cid:3) The next result analyzes the case n = 3. Note that (19) implies that e goes to one faster than e . Figure 7.
Steady-state densities e i as a function of i for a EFRBM with n = 6, λ i = 1, i = 0 , . . . ,
6, for three values of r (with q = 1 /r ). Proposition 4
The steady-state densities in the EFRBM with n = 3 satisfy e ( r ) = 1 − λ λ λ ( λ + λ ) r + o( r ) ,e ( r ) = 1 − λ λ r + o( r ) ,e ( r ) = λ λ + λ + o (cid:0) r (cid:1) , (20) and R ( r ) = λ λ λ + λ r + o (cid:0) r (cid:1) . (21)Note that this implies thatlim r → e ( r ) = lim r → e ( r ) = 1 , and lim r → R ( r ) = 0 , so again as r → r goes to 0 the repelling force between existing neighbors is veryweak, and the binding force when forming new neighbors is very strong, leading to theformation of traffic jams at the beginning of the lattice. Consequently, the steady-stateflow goes to zero.We now turn to consider the opposite case, that is, r → ∞ . The case r → ∞ 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 newneighbors (small q ). As we will see below, this leads to a phenomena that may be Figure 8.
Steady-state densities e i as a function of i for a EFRBM with n = 6, λ = 1 , λ = 1 . , λ = 0 . , λ = 0 . , λ = 1 . , λ = 0 .
75, and λ = 1 .
15, for threevalues of r , and q = 1 /r . The steady-state values for r = 1 ,
000 and r = 10 ,
000 cannotbe distinguished.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 λ = (cid:2) . . .
95 1 . .
75 1 . (cid:3) (cid:48) . For r = 1 (recall that q = 1 /r ), e = (cid:2) . . . . . . (cid:3) (cid:48) , R = 0 . . For r = 1 , e = (cid:2) . . . . . . (cid:3) (cid:48) , R = 0 . . For r = 10 , e = (cid:2) . . . . . . (cid:3) (cid:48) , R = 0 . . Fig. 8 depicts these steady-state values for the three r values. Note that the steady-statevalues for r = 1 ,
000 and r = 10 ,
000 cannot be distinguished. It may be observed thatthe values e j ( r ), j = 2 , ,
6, decrease to zero as r increases. In other words, in every pairof consecutive sites one density is very small. This “separation of densities” representsthe opposite of a traffic jam. This leads to a substantial increase in the output rate R as r increases. (cid:3) We now rigorously analyze the case r → ∞ for the EFRBM with n = 2 and n = 3. Example 9
Consider the EFRBM with n = 2. Expanding e in (18) as a Taylor seriesin q = 1 /r yields e = λ λ q + o ( q ) , e = λ λ + λ + o ( q ) , o lim r →∞ e = 0 , lim r →∞ e = λ λ + λ , lim r →∞ R = λ λ λ + λ . Thus, in this case the density at site 2 goes to zero, and this yields a positive steady-stateoutput rate. (cid:3)
Proposition 5
The steady-state densities in the EFRBM with n = 3 satisfy e ( q ) = a + b q + o (cid:0) q (cid:1) ,e ( q ) = λ λ q + o (cid:0) q (cid:1) ,e ( q ) = a + b q + o (cid:0) q (cid:1) , (22) with a , a ∈ (0 , , and R ( q ) = λ (1 − a ) + λ (cid:18) ( a − λ λ − b (cid:19) q + o (cid:0) q (cid:1) . (23)Note that this implies thatlim r →∞ e ( r ) = 0 , and lim r →∞ R ( r ) > , so again as r → ∞ 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. In-deed, 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 dependon the presence of neighbors.To study the effect of such interactions, we introduced a new deterministic compart-mental model, the EFRBM, for the flow of particles along an ordered lattice of siteswhere the transition rates between sites depend both on properties of the lattice andon nearest-neighbor interactions between the particles. The properties of the lattice aremodeled using transition rates λ i between sites. The nearest-neighbor interactions be-tween the particles are modeled using two parameters: r that represents the tendency ofa moving particle to break from an existing neighbor, and q that represents the tendencyof 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 ofoccupations between neighboring sites. It is possible to use our framework also toderive a more complete model based on binary occupation densities and transitionsdescribed by a continuous-time master equation (see, e.g. the interesting paper [44] inwhich this was done for granular channel transport). However, in such a model thestate-variables at time t represent the probability of each configuration at time t , andthe 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. Anotherimportant advantage of the EFRBM is that it is amenable to analysis using tools fromsystems and control theory, even in the non-homogeneous case. This allows to rigorouslystudy, for example, the effect of the nearest-neighbors interactions on the steady-statebehavior of the EFRBM for any set of transition rates. Our results show that suitableforces between nearby particles can greatly increase the output rate, and reveal thatthe 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-stateensity is such that any second site is empty. This represents the “opposite” of a trafficjam, and increases the steady-state flow.The phenomenological model introduced here may prove useful for other applicationsas well. For example, an important problem in vehicular traffic is to understand howhuman drivers react to nearby cars. One may also consider implementing appropriatenearest-neighbor dynamics in algorithms that control autonomous vehicles in order toreduce traffic jams and increase the flow. Of course, implementing this with a verylarge r (or q ) means very high effective transition rates, but our results suggest thateven for r not much larger than one the increase in the flow is non-negligible. Anotherinteresting topic for further research is generalizing the EFRBM to include the possibilityof attachment/detachment of particles from intermediate sites in the lattice (see [45] forsome 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 Bioin-formatics at Tel Aviv University. The research of MM is partially supported by researchgrants from the Israeli Ministry of Science, Technology & Space, the US-Israel Bina-tional 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 followsimmediately from the equations of the EFRBM. Let η i ( t ) := λ i (1 + ( q − z i +2 ( t ))(1 + ( r − z i − ( t )) , i = 0 , . . . , n, (24)with the z i s defined in (2). By (3), the EFRBM can be written as˙ x i ( t ) = η i − ( t ) x i − ( t )(1 − x i ( t )) − η i ( t ) x i ( t )(1 − x i +1 ( t )) . (25)This is just the RFM (see (1)), but with time-varying rates η i ( t ). Let a i := min { , q } min { , r } λ i ,and b i := max { , q } max { , r } λ i . It follows from (24) that a i ≤ η i ( t ) ≤ b i for all i andfor all t ≥
0. Note that for r, q > a i is strictly positive. In other words, all thetime-varying rates are uniformly separated from zero and uniformly bounded. Now theproof of Proposition 1 follows from the results in [28]. (cid:3) Proof of Proposition 2.
Combining the representation in (25) with the uniformboundedness of the rates, Proposition 1, and the results in [43] imply that the EFRBMis contractive after a small overshoot and short transient (SOST) on C n . Also, Proposi-tion 4 in [43] implies that for the EFRBM the properties of SOST and SO are equivalent,and this completes the proof. (cid:3) Proof of Fact 2.
Consider the EFRBM with n = 2 and q = 1 /r . Then (15) becomes λ (1 − e )(1 + ( 1 r − e ) = λ e (1 − e )= λ e (1 + ( r − e ) . This yields e = λ e λ + ( λ (1 − r ) − λ ) e (26)and a e + a e + λ = 0 , ith a defined in (17) and a := − λ − λ − a . The feasible solution (i.e. the onesatisfying e , e ∈ (0 ,
1) for any set of parameter values) is given by e = − a − (cid:112) a − a λ a , and (26). (cid:3) Proof of Prop. 3.
To emphasize the dependence on the parameters, write the EFRBMas ˙ x = f ( x ; v ), where v := (cid:2) λ . . . λ n r q (cid:3) (cid:48) . Note that f is an analytic function.Then the steady-state satisfies the relation f ( e ; v ) = 0. The Jacobian matrix of thisrelation with respect to x is J ( x ; v ) := ∂∂x f ( x ; v ) , which is just the Jacobian of the dynamics. Fix v ∈ R n +3++ and let e ∈ Int( C n ) denotethe corresponding steady-state, that is, f ( e ; v ) = 0 and e = h ( v ). Suppose thatthere exists a matrix measure µ : R n × n → R such that µ ( J ( e ; v )) <
0. This impliesin particular that J ( e , v ) is Hurwitz (see e.g. [46]), so it is not singular and invokingthe implicit function theorem implies that the mapping h is analytic. It follows fromthe results in [28] that such a matrix measure µ indeed exists, and this completes theproof. (cid:3) Proof of Prop. 4.
Expand e i , i = 1 , ,
3, as e i = a i + b i r + c i r + o( r ) . (27)Recall that the steady-state equations are given by R ( e ) = g ( e ) = g ( e ) = · · · = g ( e ),with the g i s given in (4). Substituting (27) yields g ( e ) = λ (1 − a ) a r + . . . ,g ( e ) = λ (1 − a ) a r + . . . ,g ( e ) = λ ( a − a ( a −
1) + . . . ,g ( e ) = λ (1 − a ) a + . . . . Since R ( e ) is bounded, we conclude that(1 − a ) a = (1 − a ) a = 0 . (28)Assume for the moment that a = 0. Then g ( e ) = λ ( a − a − b r + . . . , (29) g ( e ) = λ a + . . . , and this implies that a = 0. Now, g ( e ) = λ (1 + b ) + . . . and combining thiswith (29) yields b = −
1. Thus, e = a + b r + c r + o( r ) = − r + o (cid:0) r (cid:1) , and thisis a contradiction as e ( r ) will be strictly negative for any r > a (cid:54) = 0, so (28) yields a = 1, and also (1 − a ) a = 0. Supposethat a = 0. Then g ( e ) = − λ a b + . . . ,g ( e ) = λ (1 + b )(1 − a ) + . . . ,g ( e ) = λ a ( b − r + . . . ,g ( e ) = λ (1 − a ) b r + . . . . t follows that a b = (1 + b )(1 − a ) = 0. Since we already know that a (cid:54) = 0, b = 0. The case b = − e ( r ) < r > a = 1. But then R ( e ) = g ( e ) = − λ r + . . . and this is a contradiction. We concludethat a (cid:54) = 0, so (28) yields a = 1. Summarizing, we have a = a = 1. Now, g ( e ) = − λ b + . . . ,g ( e ) = − λ b a + . . . ,g ( e ) = λ ( a − b − r + . . . ,g ( e ) = λ (1 − b ) a r + . . . . This gives b = 0 and b a = 0. Since we already know that a (cid:54) = 0, b = 0. Now, g ( e ) = − λ c r + . . . ,g ( e ) = − λ a c r + . . . ,g ( e ) = λ (1 − a ) r + . . . ,g ( e ) = λ a r + . . . . Equating the coefficients here yields a = λ λ + λ , c = − λ λ λ ( λ + λ ) , and c = − λ /λ .Since we know that the steady-state equations admit a unique solution this yields (20),and the equation R ( e ) = g ( e ) yields (21). (cid:3) Proof of Prop. 5.
Expand e i , i = 1 , ,
3, as e i = a i + b i q + c i q + o( q ) . (30)Recall that the steady-state equations are given by R ( e ) = g ( e ) = g ( e ) = · · · = g ( e ),with the g i s given in (4). Substituting (30) yields g ( e ) = λ ( a − a −
1) + . . . ,g ( e ) = λ a ( a − a −
1) + . . . ,g ( e ) = λ a a (1 − a ) q + . . . ,g ( e ) = λ a a q + . . . . This implies that a a (1 − a ) = a a = 0 . (31)Assume for the moment that a (cid:54) = 0. Then a = a = 0. This yields g ( e ) = λ (1 − a ) + . . . ,g ( e ) = λ (1 − a ) b q + . . . ,g ( e ) = λ a (1 + b ) + . . . . This implies that a = 1 and b = −
1. This yields e ( r ) = − r + o (cid:0) r (cid:1) which is acontradiction.We conclude that a = 0. Now, g ( e ) = λ (1 − a ) + . . . ,g ( e ) = λ a (1 − a ) + . . . ,g ( e ) = λ a (1 − a ) b + . . . ,g ( e ) = λ a (1 + b ) + . . . . Equating the coefficients here yields the following. First, a (cid:54) = 0, and since e ( q ) = a + . . . , this implies that a >
0. Second, if a = 1 then a = 1 and b = − e ( q ) = − q + o ( q ). Thus, a (cid:54) = 1 and this implies that a (cid:54) = 1. Weconclude that a , a ∈ (0 , g and g yield b = λ /λ . Thisproves (22). Expanding g up to order one in q , and using R = g yields (23). (cid:3) eferences
1. Alberts B, Johnson A, Lewis J, Raff M, Roberts K, Walter P. Molecular Biologyof the Cell. 5th ed. New York: Garland Science; 2007.2. Kolomeisky AB. Motor Proteins and Molecular Motors. CRC Press; 2015.3. Schadschneider A, Chowdhury D, Nishinari K. Stochastic Transport in ComplexSystems: From Molecules to Vehicles. Elsevier; 2011.4. Kolomeisky AB. Motor proteins and molecular motors: how to operate machinesat the nanoscale. Journal of Physics: Condensed Matter. 2013;25(46):463101.5. Chowdhury D, Schadschneider A, Nishinari K. Physics of transport and trafficphenomena in biology: from molecular motors and cells to organisms. Physics ofLife Reviews. 2005; p. 318–352.6. Mallik R, Rai AK, Barak P, Rai A, Kunwar A. Teamwork in microtubule motors.Trends in Cell Biology. 2013;23:575–582.7. Blythe RA, Evans MR. Nonequilibrium steady states of matrix-product form: asolver’s guide. J Phys A Math Gen. 2007;40(46):R333–R441.8. Zia R, Dong J, Schmittmann B. Modeling translation in protein synthesis withTASEP: A tutorial and recent developments. J Stat Phys. 2011;144:405–428.9. Kriecherbauer T, Krug J. A pedestrian’s view on interacting particle sys-tems, KPZ universality and random matrices. J Phys A: Math Theor.2010;43(40):403001.10. Romano MC, Thiel M, Stansfield I, Grebogi C. Queuing phase transition: Theoryof translation. Phys Rev Lett. 2009;102:198104.11. Ciandrini L, Stansfield I, Romano MC. Ribosome traffic on mRNAs maps togene ontology: Genome-wide quantification of translation initiation rates andpolysome size regulation. PLOS Computational Biology. 2013;9(1):e1002866.12. Lakatos G, Chou T. Totally asymmetric exclusion processes with particles ofarbitrary size. J Phys A: Math Gen. 2003;36:20272041.13. Shaw LB, Zia RKP, Lee KH. Totally asymmetric exclusion process with extendedobjects: a model for protein synthesis. Phys Rev E. 2003;68:021910.14. Shaw LB, Kolomeisky AB, Lee KH. Local inhomogeneity in asymmetric simpleexclusion processes with extended objects. J Phys A Math Gen. 2004;37(6):2105.15. Srinivasa S, Haenggi M. A statistical mechanics-based framework to analyze adhoc networks with random access. IEEE Trans Mobile Computing. 2012;11:618–630.16. Turci F, Parmeggiani A, Pitard E, Romano MC, Ciandrini L. Transport on alattice with dynamical defects. Phys Rev E. 2013;87(1):012705.17. Pinkoviezky I, Gov NS. Modelling interacting molecular motors with an internaldegree of freedom. New J of Physics. 2013;15(2):025009.18. Antal T, Sch¨utz G. Asymmetric exclusion process with next-nearest-neighbor in-teraction: Some comments on traffic flow and a nonequilibrium reentrance tran-sition. Phys Rev E. 2000;62(1):83.9. Hager J, Krug J, Popkov V, Sch¨utz G. Minimal current phase and universalboundary layers in driven diffusive systems. Phys Rev E. 2001;63(5):056110.20. Teimouri H, Kolomeisky AB, Mehrabiani K. Theoretical analysis of dy-namic processes for interacting molecular motors. J Phys A: Math Theor.2015;48(6):065001.21. Celis-Garza D, Teimouri H, Kolomeisky AB. Correlations and symmetry of inter-actions influence collective dynamics of molecular motors. Journal of StatisticalMechanics: Theory and Experiment. 2015;2015(4):P04013.22. Pinkoviezky I, Gov NS. Traffic jams and shocks of molecular motors inside cellularprotrusions. Phys Rev E. 2014;89:052703.23. Leduc C, Padberg-Gehle K, Varga V, Helbing D, Diez S, Howard J. Molecularcrowding creates traffic jams of kinesin motors on microtubules. Proceedings ofthe National Academy of Sciences. 2012;109(16):6100–6105.24. Mandelkow E, Mandelkow EM. Kinesin motors and disease. Trends in CellBiology. 2002;12:585–591.25. Reuveni S, Meilijson I, Kupiec M, Ruppin E, Tuller T. Genome-scale analysis oftranslation elongation with a ribosome flow model. PLOS Computational Biology.2011;7:e1002127.26. Margaliot M, Tuller T. Stability analysis of the ribosome flow model. IEEE/ACMTrans Computational Biology and Bioinformatics. 2012;9:1545–1552.27. Smith HL. Monotone Dynamical Systems: An Introduction to the Theory ofCompetitive and Cooperative Systems. vol. 41 of Mathematical Surveys andMonographs. Providence, RI: Amer. Math. Soc.; 1995.28. Margaliot M, Sontag ED, Tuller T. Entrainment to periodic initiation andtransition rates in a computational model for gene translation. PLoS ONE.2014;9(5):e96039.29. Margaliot M, Tuller T. On the steady-state distribution in the homogeneous ribo-some flow model. IEEE/ACM Trans Computational Biology and Bioinformatics.2012;9:1724–1736.30. Zarai Y, Margaliot M, Tuller T. Explicit expression for the steady-statetranslation rate in the infinite-dimensional homogeneous ribosome flow model.IEEE/ACM Trans Computational Biology and Bioinformatics. 2013;10:1322–1328.31. Margaliot M, Tuller T. Ribosome flow model with positive feedback. J RoyalSociety Interface. 2013;10:20130267.32. Poker G, Zarai Y, Margaliot M, Tuller T. Maximizing protein translation ratein the nonhomogeneous ribosome flow model: a convex optimization approach. JRoyal Society Interface. 2014;11(100).33. Poker G, Margaliot M, Tuller T. Sensitivity of mRNA translation. Sci Rep.2015;5:12795.34. Raveh A, Zarai Y, Margaliot M, Tuller T. Ribosome flow model on a ring.IEEE/ACM Trans Computational Biology and Bioinformatics. 2015;12(6):1429–1439.5. Zarai Y, Margaliot M, Sontag ED, Tuller T. Controlling mRNA translation.2016;.36. Zarai Y, Margaliot M, Tuller T. Optimal down regulation of mRNA translation.Sci Rep. 2017;7(41243).37. Zarai Y, Margaliot M, Tuller T. On the ribosomal density thatmaximizes protein translation rate. PLOS ONE. 2016;11(11):1–26.doi:10.1371/journal.pone.0166481.38. Teimouri H, Kolomeisky AB, Mehrabiani K. Theoretical analysis of dy-namic processes for interacting molecular motors. J Phys A: Math Theor.2015;48(6):065001.39. Nadler W, Schulten K. Generalized moment expansion for observables of stochas-tic processes in dimensions d >d >