A Deterministic Mathematical Model for Bidirectional Excluded Flow with Langmuir Kinetics
aa r X i v : . [ q - b i o . S C ] S e p A Deterministic Mathematical Model forBidirectional Excluded Flow withLangmuir Kinetics
Yoram Zarai, Michael Margaliot, and Tamir Tuller
Abstract
In many important cellular processes, including mRNA translation, gene transcription, phosphotransfer, andintracellular transport, biological “particles” move along some kind of “tracks”. The motion of these particles canbe modeled as a one-dimensional movement along an ordered sequence of sites. The biological particles (e.g.,ribosomes, RNAPs, phosphate groups, motor proteins) have volume and cannot surpass one another. In some cases,there is a preferred direction of movement along the track, but in general the movement may be two-directional,and furthermore the particles may attach or detach from various regions along the tracks (e.g. ribosomes may dropoff the mRNA molecule before reaching a stop codon).We derive a new deterministic mathematical model for such transport phenomena that may be interpreted as thedynamic mean-field approximation of an important model from mechanical statistics called the asymmetric simpleexclusion process (ASEP) with Langmuir kinetics. Using tools from the theory of monotone dynamical systemsand contraction theory we show that the model admits a unique equilibrium, and that every solution converges tothis equilibrium. This means that the occupancy in all the sites along the lattice converges to a steady-state valuethat depends on the parameters but not on the initial conditions. Furthermore, we show that the model entrains (orphase locks) to periodic excitations in any of its forward, backward, attachment, or detachment rates.We demonstrate an application of this phenomenological transport model for analyzing the effect of ribosomedrop off in mRNA translation. One may perhaps expect that drop off from a jammed site may increase the totalflow by reducing congestion. Our results show that this is not true. Drop off has a substantial effect on the flow,yet always leads to a reduction in the steady-state protein production rate.
Index Terms
Monotone dynamical systems, systems biology, synthetic biology, mRNA translation, gene transcription, ri-bosome flow model, ribosome drop off, Langmuir kinetics, bi-directional flow, intracellular transport, contraction
The research of MM and TT is partially supported by research grants from the Israeli Ministry of Science, Technology & Space and theBinational Science Foundation. The research of MM is also supported by a research grant form the Israeli Science Foundation.Y. Zarai is with the School of Electrical Engineering, Tel-Aviv University, Tel-Aviv 69978, Israel. E-mail: [email protected]. Margaliot is with the School of Electrical Engineering and the Sagol School of Neuroscience, Tel-Aviv University, Tel-Aviv 69978,Israel. E-mail: [email protected]. Tuller is with the department of Biomedical Engineering and the Sagol School of Neuroscience, Tel-Aviv University, Tel-Aviv 69978,Israel. E-mail: [email protected] theory, contraction after a short transient, entrainment.
I. I
NTRODUCTION
Movement is essential for the functioning of cells. Cargoes like organelles and vesicles must be carriedbetween different locations in the cells. The information encoded in DNA and mRNA molecules must bedecoded by “biological machines” (RNA polymerases and ribosomes) that move along these moleculesin a sequential order.Many of these important biological transport processes are modeled as the movement of particles alongan ordered chain of sites. In the context of intercellular transport, the particles are motor proteins andthe chain models actin filaments or microtubules. In transcription, the particles are RNAPs moving alongthe DNA molecule, and in translation the particles are ribosomes moving along the mRNA molecule (seeFigure 1).The movement in such processes may be unidirectional, as in mRNA translation elongation, or bidi-rectional, as in transcription or translation initiation. Indeed, the normal forward flow of the RNAP maybe interrupted, due to transcription errors and various obstacles such as nucleosomes, in which case theRNAP tracks back a few nucleotides and then resumes its normal forward flow [54], [41], [8], [10].Translation initiation in eukaryotes usually includes diffusion from the 5’end of the transcript towardsthe start codon [1]. This diffusion process is believed to be bi-directional, but with a preference to the5’ →
3’ direction. The movement of motor proteins like kinesin and dynein along microtubules is typicallyunidirectional, but can be two-directional as well [1].To increase efficiency, many particles may move simultaneously along the same track thus pipeliningthe production process. For example, to increase translation efficiency, a number of ribosomes may actsimultaneously as polymerases on the same mRNA molecule [66], [4].The moving biological particles have volume and usually cannot overtake a particle in front of them.This means that a slowly moving particle may lead to the formation of a traffic jam behind it. For example,Leduc et al. [31] have studied Kip3, a yeast kinesin-8 family motor, and demonstrated that motor proteintraffic jams can exist, given the right conditions. Other studies have suggested that traffic jams of RNAP[ribosomes] may evolve during transcription [translation] [4], [27], [9].In some of these biological transport processes, the biological machines may either attach or detachat various sites along the tracks. For example, ribosomes may detach from the mRNA molecule before
A.B.C.
Initiation Elongation
TranscriptionTranslationTransport
Microtubule Slow codon Nucleosome
Fig. 1. Biological processes that can be studied using the model suggested in this paper. reaching the stop codon due to traffic “jams” and ribosome-ribosome interactions or due to depletion inthe concentration of tRNAs [72], [28], [57]. Also, it is known that kinesin-family motor proteins are moresusceptible to dissociation when their pathway is blocked [14], [62].Defects in these transport processes may lead to severe diseases or may even be lethal. For example, [53]lists the implications of malfunctions of protein motors in disease and developmental defects.Developing a better understanding of these dynamical biological processes by combining mathematicalmodeling and biological experiments will have far reaching implications to basic science in fields such asmolecular evolution and functional genomics, as well as applications in synthetic biology, biotechnology, human health, and more. Mathematical or computational modeling is especially important in developingapproaches for manipulating and controlling these processes, e.g. in order to optimize various goals inbiotechnology.A standard model for such transport processes is the asymmetric simple exclusion process (ASEP) [55],[73]. This is a stochastic model describing particles that hop along an ordered lattice of sites. Each sitecan be either empty or occupied by a single particle, and a particle can only hop to an empty site. This“simple exclusion principle” represents the fact that the particles have volume and cannot overtake oneanother. Simple exclusion generates an indirect coupling between the particles. In particular, traffic jamsmay develop behind a slow-moving particle.In ASEP, a particle may hop to any of the two neighboring sites (but only if they are free). Typically,a particle can attach the lattice in one of its ends and detach from the other end. When particles canalso attach or detach at internal sites along the lattice, the model is referred to as ASEP with
Langmuirkinetics . In the special case where the hops are unidirectional, ASEP is sometimes referred to as the totallyasymmetric simple exclusion process (TASEP). A TASEP-like system with Langmuir kinetics has beenused to model limit order markets in [65], and is often used in modeling molecular motor traffic [42], [43],[32], [33], [16]. More generally, ASEP has become a fundamental model in non-equilibrium statisticalmechanics, and has been applied to model numerous natural and artificial processes including traffic andpedestrian flow, the movement of ants, evacuation dynamics, and more [52].In this paper, we introduce a deterministic mathematical model that may be interpreted as the dynamic mean-field approximation of ASEP with Langmuir kinetics (MFALK). We analyze the MFALK using toolsfrom systems and control theory. In particular, we apply some recent developments in contraction theoryto prove that the model is globally asymptotically stable, and that it entrains to periodic excitations inthe transition/attachment/detachment rates. In other words, if these rates change periodically in time withsome common period T then all the state-variables in the MFALK converge to a periodic solution withperiod T . This is important because many biological processes are excited by periodic signals (e.g. the24h solar day or the periodic cell-division process), and proper functioning requires phase-locking orentrainment to these excitations.Our work is motivated by the analysis of a model for mRNA translation called the ribosome flowmodel (RFM) [48]. This is the mean-field approximation of the unidirectional TASEP without
Langmuirkinetics (see, e.g., [52, section 4.9.7] and [6, p. R345]). Recently, the RFM has been studied extensively using tools from systems and control theory [36], [71], [37], [38], [35], [44], [45], [47], [70]. The analysisis motivated by implications to many important biological questions. For example, the sensitivity of theprotein production rate to the initiation and elongation rates along the mRNA molecule [45], maximizationof protein production rate [44], the effect of ribosome recycling [38], [47], and the consequences ofcompetition for ribosomes on large-scale simultaneous mRNA translation in the cell [46] (see also [19],[2] for some related models).The MFALK presented here is much more general than the RFM, and can thus be used to model andanalyze many transport phenomena, including all the biological processes mentioned above, that cannotbe captured using the RFM. We demonstrate this by using the MFALK to model and analyze mRNAtranslation with ribosome drop off - a feature that cannot be modeled using the RFM.Ribosome drop off is a fundamental phenomena that has received considerable attention (see, e.g., [57],[25], [24], [68], [7], [61], [18], [28], [22], [20]). In many cases, ribosome drop off is deleterious to thecell since translation is the most energetically consuming process in the cell and, furthermore, drop offyields truncated, non-functional proteins. Thus, transcripts undergo selection to minimize drop off or itsenergetic cost [67], [63], [61], [18], [28], [20]. There are various hypotheses on the biological advantagesof ribosome drop off. For example, Zaher and Green [69] have suggested that ribosome drop off is relatedto proof reading. One may perhaps expect that another advantage is that drop off from a jammed sitemay increase the total flow by reducing congestion. Our results using analysis of the MFALK show thatthis is not true. Drop off has a substantial effect on the flow, yet it always leads to a reduction in thesteady-state protein production rate.The remainder of this paper is organized as follows. The next section describes the new mathematicalmodel. Section III presents our main analysis results. Section IV describes the application of the MFALKto model mRNA translation with ribosome drop off. The final section concludes and describes possibledirections for further research. To streamline the presentation, all the proofs are placed in the Appendix.II. T
HE MODEL
The MFALK is a set of n first-order nonlinear differential equations, where n denotes the numberof compartments or sites along the “track”. Each site is associated with a state variable x i ( t ) ∈ [0 , describing the normalized “level of occupancy” at site i at time t , with x i ( t ) = 0 [ x i ( t ) = 1 ] representingthat site i is completely free [full] at time t . Since x i ( t ) ∈ [0 , for all t , it may also be interpreted asthe probability that site i is occupied at time t . The MFALK contains four sets of non-negative parameters: • λ i , i = 0 , . . . , n , controls the forward transition rate from site i to site i + 1 , • γ i , i = 0 , . . . , n , controls the backward transition rate from site i + 1 to site i , • β i , i = 1 , . . . , n , controls the attachment rate to site i , • α i , i = 1 , . . . , n , controls the detachment rate from site i ,where we arbitrarily refer to left-to-right flow along the chain as forward flow, and to flow in the otherdirection as backward flow.The dynamical equations describing the MFALK are: ˙ x = λ (1 − x ) + γ x (1 − x ) + β (1 − x ) − λ x (1 − x ) − γ x − α x , ˙ x = λ x (1 − x ) + γ x (1 − x ) + β (1 − x ) − λ x (1 − x ) − γ x (1 − x ) − α x , ... ˙ x n − = λ n − x n − (1 − x n − ) + γ n − x n (1 − x n − ) + β n − (1 − x n − ) − λ n − x n − (1 − x n ) − γ n − x n − (1 − x n − ) − α n − x n − , ˙ x n = λ n − x n − (1 − x n ) + γ n (1 − x n ) + β n (1 − x n ) − λ n x n − γ n − x n (1 − x n − ) − α n x n . (1)To explain these equations, consider for example the equation for the change in the occupancy in site ,namely, ˙ x = λ x (1 − x ) + γ x (1 − x ) + β (1 − x ) − λ x (1 − x ) − γ x (1 − x ) − α x . The term λ x (1 − x ) represents the flow from site to site . This increases with the occupancy insite , and decreases with the occupancy in site . In particular, this term becomes zero when x = 1 ,i.e. when site is completely full. This is a “soft” version of the hard exclusion principle in ASEP: theeffective entry rate into a site decreases as it becomes fuller. Note that the constant λ ≥ describes themaximal possible transition rate from site to site . Similarly, the term λ x (1 − x ) represents the flowfrom site to site . The term γ x (1 − x ) [ γ x (1 − x ) ] represents the backward flow from site tosite [site to site ]. Note that these terms also model soft exclusion. The term β (1 − x ) representsattachment of particles from the environment to site , whereas α x represents detachment of particlesfrom site to the environment (see Fig. 2).The MFALK is a compartmental model [21], [51], as every state-variable describes the occupancy in a Output rate
Fig. 2. Topology of the MFALK. compartment (e.g., a site along the the mRNA, gene, microtubule), and the dynamical equations describethe flow between these compartments and the environment. Compartmental models play an importantrole in pharmacokinetics, enzyme kinetics, basic nutritional processes, cellular growth, and pathologicalprocesses, such as tumourigenesis and atherosclerosis (see, e.g., [21], [17] and the references therein). Morespecifically, the MFALK is a nonlinear tridiagonal compartmental model, as every ˙ x i directly dependson x i − , x i , and x i +1 only.Note also that n X i =1 ˙ x i = λ (1 − x ) − γ x + β (1 − x ) − α x + γ n (1 − x n ) − λ n x n + β n (1 − x n ) − α n x n + n − X i =2 ( β i (1 − x i ) − α i x i ) . (2)The term on the right-hand side of the first [second] line here represents the change in x [ x n ] due to theflow between the environment and site [site n ], whereas the term on the third line represents the flowbetween internal sites and the environment.The output rate from site n at time t is the total flow from this site to the environment: R ( t ) : = ( λ n + α n ) x n ( t ) − ( γ n + β n )(1 − x n ( t )) . (3)Note that R ( t ) may be positive, zero, or negative.In the particular case where α i = β i = γ i = 0 for all i the MFALK becomes the RFM, i.e. thedynamic mean-field approximation of the unidirectional TASEP with open boundary conditions andwithout Langmuir kinetics. Let x ( t, a ) denote the solution of (1) at time t ≥ for the initial condition x (0) = a . Since thestate-variables correspond to normalized occupancy levels, we always assume that a belongs to the closed n -dimensional unit cube: C n := { x ∈ R n : x i ∈ [0 , , i = 1 , . . . , n } . Let int( C n ) denote the interior of C n , and let ∂C n denote the boundary of C n . The next section analyzesthe MFALK defined in (1). III. M AIN R ESULTS
A. Invariance and persistence
It is straightforward to show that C n is an invariant set for the dynamics of the MFALK, that is,if a ∈ C n then x ( t, a ) ∈ C n for all t ≥ . The following result shows that a stronger property holds.Recall that all the proofs are placed in the Appendix. For notational convenience, let α := 0 , γ := 0 , α n +1 := 0 , and β n +1 := 0 . Proposition 1
Suppose that at least one of the following two conditions holds: λ i + β i +1 > , for all i ∈ { , . . . , n } , (4) or γ i + α i +1 > , for all i ∈ { , . . . , n } . (5) Then for any τ > there exists d = d ( τ ) ∈ (0 , / such that d ≤ x i ( t + τ, a ) ≤ − d, (6) for all a ∈ C n , all i ∈ { , . . . , n } , and all t ≥ . This means in particular that trajectories that emanate from the boundary of C n “immediately” enter C n .This result is useful because as we will see below on the boundary of C n the MFALK looses someimportant properties. For example, the Jacobian matrix of the dynamics (1) is irreducible on int( C n ) , butbecomes reducible on some points on the boundary of C n . B. Contraction
Differential analysis and in particular contraction theory proved to be a powerful tool for analyzingnonlinear dynamical systems. In a contractive system, trajectories that emanate from different initialconditions contract to each other at an exponential rate [34], [49], [3]. Let | · | : R n → R + denote the L norm, i.e. for z ∈ R n , | z | = | z | + · · · + | z n | . Proposition 2
Let η := max {− λ − γ − α − β , − α − β , . . . , − α n − − β n − , − λ n − γ n − α n − β n } . Note that η ≤ . For any a, b ∈ C n and any t ≥ , | x ( t, a ) − x ( t, b ) | ≤ exp( ηt ) | a − b | . (7)This implies that the L distance between any two trajectories contracts with the exponential rate η .Roughly speaking, this also means that increasing all the sums α i + β i , i = 1 , . . . , n , makes the system“more contractive”. Indeed, these parameters have a direct stabilizing effect on the dynamics of site i ,whereas the other parameters affect the site indirectly via the coupling to the two adjacent sites.When η = 0 , (7) only implies that the L distance between trajectories does not increase. This propertyis not strong enough to prove the asymptotic properties described in the subsections below. Indeed, in thiscase it is possible that the MFALK will not be contractive with respect to any fixed norm. Fortunately, acertain generalization of contraction turns out to hold in this case.Consider the time-varying dynamical system ˙ x ( t ) = f ( t, x ( t )) , (8)whose trajectories evolve on a compact and convex set Ω ⊂ R n . Let x ( t, t , a ) denote the solution of (8) attime t for the initial condition x ( t ) = a . System (8) is said to be contractive after a small overshoot (SO)[39] on Ω w.r.t. a norm | · | : R n → R + if for any ε > there exists ℓ = ℓ ( ε ) > such that | x ( t, t , a ) − x ( t, t , b ) | ≤ (1 + ε ) exp( − ℓt ) | a − b | , for all a, b ∈ Ω and all t ≥ t ≥ . Intuitively speaking, this means contraction with an exponential rate,but with an arbitrarily small overshoot of ε . Proposition 3
Suppose that λ i + γ i > , for all i ∈ { , . . . , n − } , (9) and that at least one of the two conditions (4) , (5) holds. Then the MFALK is SO on C n w.r.t. the L norm, that is, for any ε > there exists ℓ = ℓ ( ε ) > such that | x ( t, a ) − x ( t, b ) | ≤ (1 + ε ) exp( − ℓt ) | a − b | , (10) for all a, b ∈ C n and all t ≥ . Note that if λ i + γ i = 0 for some i ∈ { , . . . , n − } , that is λ i = γ i = 0 , then the MFALK decouplesinto two separate MFALKs: one containing sites , . . . , i , and the other containing sites i + 1 , . . . , n . Thus,assuming (9) incurs no loss of generality.There is an important difference between Propositions 2 and 3. If η < then Proposition 2 providesan explicit exponential contraction rate. If η = 0 then Proposition 3 can be used to deduce SO, but in thisresult the contraction rate ℓ depends on ε and is not given explicitly.The contraction results above imply that the MFALK satisfies several important asymptotic properties.These are described in the following subsections. C. Global asymptotic stability
Since the compact and convex set C n is an invariant set of the dynamics, it contains an equilibriumpoint e . By Proposition 1, e ∈ int( C n ) . Applying (10) with b = e yields the following result. Corollary 1
Suppose that the conditions in Proposition 3 hold. Then the MFALK admits a uniqueequilibrium point 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 the rates determine a unique distribution profile along the lattice, and that all trajectoriesemanating from different initial conditions in C n asymptotically converge to this distribution. In addition,perturbations in the occupancy levels along the sites will not change this asymptotic behavior of thedynamics. This also means that various numerical solvers of ODEs will work well for the MFALK (seee.g. [13]). Example 1
Fig. 3 depicts the trajectories of a MFALK with n = 3 , λ = 1 . , λ = 1 . , λ = 0 . , λ = 0 . , γ i = λ i − . , i = 0 , . . . , , α = 0 , α = 0 . , α = 0 , β = 0 , β = 0 . , β = 0 , for six initialconditions in C n . It may be seen that all trajectories converge to an equilibrium point e ∈ int( C ) . (cid:3) x x x Fig. 3. Trajectories of the MFALK in Example 1 for six initial conditions in C . The MFALK (1) can be written as ˙ x i = f i − ( x ) − f i ( x ) + g i ( x i ) , i = 1 , . . . , n, (11)where f ( x ) := λ (1 − x ) − γ x ,f i ( x ) := λ i x i (1 − x i +1 ) − γ i x i +1 (1 − x i ) , i = 1 , . . . , n − ,f n ( x ) := λ n x n − γ n (1 − x n ) ,g i ( x i ) := β i (1 − x i ) − α i x i , i = 1 , . . . , n. (12)At steady-state, i.e. for x = e , the left-hand side of all the equations in (11) is zero, so f i − ( e ) = f i ( e ) − g i ( e i ) , i = 1 , . . . , n. (13)Let v := h α , . . . , α n , β , . . . , β n , γ , . . . , γ n , λ , . . . , λ n i ′ ∈ R n +2+ denote the parameters of the MFALK.It follows from (13) that if we multiply all these parameters by c > then e will not change, that is, e ( cv ) = e ( v ) . Let R := ( λ n + α n ) e n − ( γ n + β n )(1 − e n ) , (14)denote the steady-state output rate . Then R ( cv ) = cR ( v ) , for all c > , that is, the steady-state productionrate is homogeneous of order one w.r.t. the parameters. By (13), R = f n ( e ) − g n ( e n )= f i ( e ) + n − X j = i +1 g j ( e j ) , i = 0 , . . . , n − . (15)This yields the following set of recursive equations relating the steady-state occupancy levels and theoutput rate in the MFALK: e n = R + γ n + β n λ n + γ n + β n + α n ,e i = R + γ i e i +1 − P n − j = i +1 g j ( e j ) λ i (1 − e i +1 ) + γ i e i +1 , i = n − , . . . , , (16)and also e = λ + β − R + P n − j =2 g j ( e j ) λ + γ + β + α . For a given v , this is a set of n + 1 equations in the n + 1 unknowns: e , . . . , e n , R . Example 2
Consider the MFALK with dimension n = 2 . Then (16) becomes e = R + γ + β λ + γ + α + β ,e = R + γ e λ (1 − e ) + γ e , (17)and also e = λ + β − Rλ + γ + β + α . This yields the polynomial equation a R + a R + a = 0 , where a := λ − γ ,a := ( λ − γ )( γ + β − λ − β ) − λ z − z z − z γ ,a := ( λ + β ) λ ( λ + α ) − ( γ + α ) γ ( γ + β ) , with z := λ + γ + α + β and z := λ + γ + α + β .Note that the polynomial equation admits several solutions R , but only one solution corresponds tothe unique equilibrium point e ∈ C . For example, for λ i = 1 , γ i = 2 , β i = 3 , and α i = 4 for all i thepolynomial equation becomes − R − R −
40 = 0 . This admits two solutions R = ( − s − / and R = (3 s − / , with s := √ . Substituting R in (17) yields e = [ e e ] ′ , with e < , sothis is not a feasible solution. Substituting R in (17) yields (all numbers are to four digit accuracy) e = h . . i ′ ∈ C , which is the unique feasible solution. Thus, the steady-state output rate is R = − . . (cid:3) In general, (16) can be transformed into a polynomial equation for R . The next result shows that thedegree of this polynomial equation grows quickly with n . Proposition 4
Consider the MFALK with dimension n and with λ i = γ i , α i = 0 , β i = 0 , for all i .Then generically Eq. (16) may be written as w ( R ) = 0 , where w ( R ) is a polynomial equation in R ofdegree ⌊ n ⌋ , and with coefficients that are algebraic functions of the rates. We note that this is exponential increase in the degree of the polynomial equation is a feature of theMFALK that does not take place in the RFM. Indeed, in the RFM the degree of the polynomial equationfor the steady-state production rate grows linearly with n .Let sgn( · ) : R → {− , , } denote the sign function, i.e. sgn( y ) = , y > , , y = 0 , − , y < . An interesting question is what is sgn( R ) . Indeed, if this is positive (negative) then this means that thereis a net steady-state flow from left to right (right to left). The next subsection describes a special casewhere this question can be answered rigorously.
1) Bidirectional flow with no Langmuir kinetics:
In the case where β i = α i = 0 , i = 1 , . . . , n , i.e. asystem with no internal attachments and detachments, Eq. (15) becomes R = f i ( e ) , i = 0 , . . . , n. (18) Proposition 5
Consider the case where α i = β i = 0 , i = 1 , . . . , n , and suppose that (9) holds. Then sgn( R ) = sgn n Y i =0 λ i − n Y i =0 γ i ! . (19) In particular, if Q ni =0 λ i = Q ni =0 γ i then R = 0 , and e i = Q i − j =0 λ j Q i − j =0 λ j + Q i − j =0 γ j = Q nj = i γ j Q nj = i γ j + Q nj = i λ j , i = 1 , . . . , n. (20)Eq. (19) means that in the case of no Langmuir kinetics the steady-state output from the right hand-sideof the chain will be positive [negative] if the the product of the forward rates is larger [smaller] thanthe product of the backward rates. In transcription and translation the steady state flow from the righthand-side of the chain should always be positive, but in other cases, e.g. transport along microtubules,the steady state flow may be either positive or negative. D. Entrainment
Assume now that some or all of the rates are time-varying periodic functions with the same period T .This may be interpreted as a periodic excitation of the system. Many biological processes are affected bysuch excitations due for example to the periodic 24h solar day or the periodic cell-cycle division process.For example, translation elongation factors, tRNAs, translation and transcription initiation factors, ATPlevels, and more may change in a periodic manner and affect various rates that appear in the MFALK.A natural question is will the state-variables of the MFALK converge to a periodic pattern with period T ?We will show that this is indeed so, i.e. the MFALK entrains to a periodic excitation in the rates. In order tounderstand what this means, consider a different setting, namely, using the MFALK to model traffic flow.Then the rates may correspond to traffic lights, changing in a periodic manner, and the state-variables arethe density of the moving particles (cars) along different sections of the road, so entrainment correspondsto what is known as the “green wave” (see e.g. [26] and the references therein).We say that a function f is T -periodic if f ( t + T ) = f ( t ) for all t . Assume that the λ i s, γ i s, α i s and β i sare uniformly bounded, non-negative, time-varying functions satisfying: • there exists a (minimal) T > such that all the λ i ( t ) s, γ i ( t ) s, α i ( t ) s, and β i ( t ) s are T -periodic. • there exist c , c > such that at least one of the following two conditions holds for all time tλ i ( t ) + β i +1 ( t ) > c , i = 0 , . . . , n, (21) γ i ( t ) + α i +1 ( t ) > c , i = 0 , . . . , n. (22) • there exists c > such that λ i ( t ) + γ i +1 ( t ) > c , i = 0 , . . . , n. (23)We refer to this model as the Periodic MFALK (PMFALK) . Theorem 1
Consider the PMFALK with dimension n . There exists a unique function φ ( · ) : R + → int( C n ) ,that is T -periodic, and for any a ∈ C n the trajectory x ( t, a ) converges to φ as t → ∞ . Thus, the PMFALK entrains (or phase-locks) to the periodic excitation in the parameters. In particular,this means that the output rate R ( t ) in (3) converges to the unique T -periodic function: ( λ n ( t ) + γ n ( t ) + β n ( t ) + α n ( t )) φ n ( t ) − γ n ( t ) − β n ( t ) . Note that since a constant function is a periodic function for all T ≥ , Theorem 1 implies that entrainmentholds also in the particular case where a single parameter is oscillating (with period T > ), while allother parameters are constant. Note also that Corollary 1 follows from Theorem 1. Example 3
Consider the MFALK with dimension n = 3 , parameters: λ ( t ) ≡ . , λ ( t ) ≡ . , λ ( t ) =1 + 0 . πt/ , λ ( t ) ≡ . , γ ( t ) ≡ . , γ ( t ) = 0 . πt/
4) + 1 / , γ ( t ) ≡ . , γ ( t ) ≡ . , α ( t ) ≡ , α ( t ) ≡ . , α ( t ) ≡ , β ( t ) ≡ , β ( t ) = 0 . πt/
2) + 1 / , β ( t ) ≡ , and initialcondition x (0) = h . . . i ′ . Note that all the rates here are periodic, with a minimal commonperiod T = 8 . Fig. 4 depicts x i ( t ) , i = 1 , , , as a function of t . It may be seen that each state variableconverges to a periodic function with period T = 8 . (cid:3) E. Strong Monotonicity
Recall that a proper cone K ⊆ R n defines a partial ordering in R n as follows. For two vectors a, b ∈ R n ,we write a ≤ b if ( b − a ) ∈ K ; a < b if a ≤ b and a = b ; and a ≪ b if ( b − a ) ∈ int( K ) . Thesystem ˙ y = f ( y ) is called monotone if a ≤ b implies that y ( t, a ) ≤ y ( t, b ) for all t ≥ . In other t Fig. 4. State variables x ( t ) [solid line]; x ( t ) [dashed line]; and x ( t ) [dotted line] as a function of t in Example 3. Note that each statevariable converges to a periodic function with a period T = 8 . words, the flow preserves the partial ordering [60]. It is called strongly monotone if a < b impliesthat y ( t, a ) ≪ y ( t, b ) for all t > .From here on we consider the particular case where the cone is K := R n + . Then a ≤ b if a i ≤ b i forall i , and a ≪ b if a i < b i for all i . A system that is monotone with respect to this partial ordering iscalled cooperative . Proposition 6
For any a, b ∈ C n , with a ≤ b , the solutions of the MFALK satisfy x ( t, a ) ≤ x ( t, b ) , for all t ≥ . (24) Furthermore, if (9) holds then x ( t, a ) ≪ x ( t, b ) , for all t > . (25)To explain this, consider two initial densities a and b with a i ≤ b i for all i , that is, b corresponds to alarger or equal density at each site. Then the trajectories x ( t, a ) and x ( t, b ) emanating from these initialconditions continue to satisfy the same relationship between the densities, namely, x i ( t, a ) ≤ x i ( t, b ) , forall i and for all time t ≥ .The MFALK is thus a strongly cooperative tridiagonal system (SCTS) on int( C n ) . Some of theproperties deduced above using contraction theory can also be deduced using this property [59]. Remark 1
Suppose that we augment the MFALK into a model of n + 1 ODEs in n + 1 state-variablesby adding to it the equation ˙ x n +1 = − λ (1 − x ) − γ x − β (1 − x ) + α x − γ n (1 − x n ) + λ n x n − β n (1 − x n ) + α n x n − n − X i =2 ( β i (1 − x i ) − α i x i ) . that is, ˙ x n +1 = − P ni =1 ˙ x i (see (2) ). Let ˜ x denote the vector of the n + 1 state-variables. Clearly, thisaugmented model admits a first integral H (˜ x ( t )) := P n +1 i =1 ˜ x i ( t ) . Also, for any initial condition in ˜ x (0) ∈ C n × R + all the state-variables remain bounded, as the first n state-variables remain in C n and ˜ x n +1 ( t ) = H (˜ x (0)) − P ni =1 ˜ x i ( t ) for all t ≥ . It is straightforward to verify that the augmented system is acooperative system, and that if (9) holds then it is a SCTS. SCTS systems that admit a non-trivial firstintegral have many desirable properties (see, e.g. [40]).F. Effect of attachment and detachment One may perhaps expect that detachment from a jammed site may increase the total flow by reducingcongestion. The next result shows that this is not so. Detachment always decreases the steady-stateproduction rate R . Similarly, attachment always increases R . Proposition 7
Consider a MFALK with dimension n . Suppose that the conditions in Proposition 3 hold.Then ∂e i ∂α j < , and ∂e i ∂β j > , for all i, j . Also, ∂R∂α j < , and ∂R∂β j > for all j = 0 , , . . . , n − . This means that an increase in any of the detachment [attachment] rates decreases [increases] thesteady-state density in all the sites. Also, an increase in any of the internal detachment [attachment] ratesdecreases [increases] the steady-state production rate. The next example demonstrates this.
Example 4
Consider the MFALK with n = 3 , λ i = 1 , γ i = 0 , i = 0 , , , , β i = α = 0 , i = 1 , , .Fig. 5 depicts R as a function of α ∈ [0 , and α ∈ [0 , . It may be seen that R decreases with both α and α . (cid:3) We note that the analytical results in Proposition 7 agree well with the simulation results obtainedusing a TASEP model for translation that included alternative initiation along the mRNA and ribosomedrop-off [74].The next section describes an application of the MFALK to a biological process. α α Fig. 5. R as a function of α ∈ [0 , and α ∈ [0 , for the MFALK in Example 4. IV. A
N APPLICATION : MODELING M
RNA
TRANSLATION WITH RIBOSOME DROP OFF
It is believed that during mRNA translation ribosome movement is unidirectional from the 5’ end tothe 3’ end, and that ribosomes do not enter in the middle of the coding regions. However, ribosomes candetach from various sites along the mRNA molecule due for example to collisions between ribosomes.This is known as ribosome drop off.As mentioned in the introduction, ribosome drop off has been the topic of numerous studies [57], [25],[24], [68], [7], [61], [18], [28], [22], [20], [29]. It was suggested that in some cases ribosome drop offis important for proof reading [69], and also that ribosome stalling/abortion plays a role in translationalregulation (e.g. see [56], [74]).It is clear that ribosome abortion has drawbacks. Indeed, translation is the most energetically consumingprocess in the cell, and abortion results in truncated, non-functional and possibly deleterious proteins. Itis believed that transcripts undergo evolutionary selection to minimize abortion and/or its energetic cost[67], [63], [61], [18], [28], [20]. Nevertheless, there seems to be a certain minimal abortion rate even innon-stressed conditions [57], [29]. This basal value was estimated (see more details below) to be of theorder or − − − abortion events per codon in E. coli . In other words, in every codon one out of , − , decoding ribosomes aborts. This value is non-negligible. If we consider a drop-off rate of ∗ − per codon along a coding region of codons (approximately the average coding region lengthfor E. coli ) then on average, around out of every ribosomes will fail to complete the translation of the mRNA.To model translation with ribosome drop off, we use the MFALK with γ i = 0 (i.e. no backwardsmotion) and β i = 0 (i.e. no attachment to internal sites along the chain) for all i . Changing the values ofthe α i s allows to model and analyze the effect of ribosome drop off at different sites along the mRNAmolecule. We assume that λ i > , for all i, (26)as otherwise the chain decouples into two smaller, disconnected chains. Note that (26) implies that theconditions in Proposition 3 hold, so the model is SO on C n w.r.t. the L norm, and thus admits a uniqueglobally asymptotically stable equilibrium point e ∈ int( C n ) .We study the effect of ribosome drop off on the steady-state protein production rate and ribosomedensity using real biological data. To this end, we considered S. cerevisiae genes (see Figures 6 and 7)with various mRNA levels (all genes were sorted according to their mRNA levels and genes wereuniformly sampled from the list). Similarly to the approach used in [48], we divided the mRNAs relatedto these genes to non-overlapping pieces. The first piece includes the first codons that are related tovarious stages of initiation [63]. The other pieces include non-overlapping codons each, except for thelast one that includes between and codons.To model the translation dynamics in these mRNAs using MFALK, we model every piece of mRNAas a site. We estimated the elongation rates λ i at each site using ribo-seq data for the codon decodingrates [12], normalized so that the median elongation rate of all S. cerevisiae mRNAs becomes . codonsper second [23]. The site rate is ( site time ) − , where site time is the sum over the decoding times ofall the codons in the piece of mRNA corresponding to this site. These rates thus depend on variousfactors including availability of tRNA molecules, amino acids, Aminoacyl tRNA synthetase activity andconcentration, and local mRNA folding [12], [1], [63].The initiation rate λ (that corresponds to the first piece) was estimated based on the ribosome densityper mRNA levels, as this value is expected to be approximately proportional to the initiation rate wheninitiation is rate limiting [48], [36]. Again we applied a normalization that brings the median initiationrate of all S. cerevisiae mRNAs to be . [9].We analyzed the effect of uniform ribosome drop off with a rate in the range of − to − per codon.This corresponds to α = · · · = α n := α c , i.e., all the α i s are equal, and α c denote their common value.Since we assumed codons per site, α c values range from − to − (ten times the rate associated α c YGR192CYOR248WYGR165WYGL223CYBR200WYEL077CYMR031CYNL095CYJR053WYER106W
Fig. 6. Reduction in steady-state mean density ρ in percent as a function of α c ∈ [10 − , − ] for S. cerevisiae genes. with a single codon). This makes sense as in the MFALK the level of occupancy in a site is related tothe probability to see a ribosome in this site.Let ρ := P ni =1 e i n , denote the steady-state mean ribosomal density. Figures 6 and 7 depict ρ and R in our model as a functionof α c ∈ [10 − , − ] . In these figures the genes in the legends are sorted according to their expressionlevels: the gene at the top (YGR192C) has the highest mRNA levels while the gene at the bottom(YER106W) has the lowest levels. It may be seen that as the drop off (detachment) rate α c increases from − to − , ρ decreases by about , and R decreases by about . This demonstrate the significantramifications that ribosomal drop off is expected to have on translation and the importance of modelingdrop off.Note also that there is a strong variability in the effect of drop off on the different genes: for mRNAswith higher expression levels (i.e. mRNAs with higher copy number in the cell) the drop off effect isweaker. It is possible that this is related to stronger evolutionary selection for lower drop off rate in geneswith higher mRNA levels. Indeed, highly expressed genes “consume” more ribosomes (due to highermRNA levels), so a given (per-mRNA) drop off rate is expected to be more deleterious to the cell, anda mutation which decreases the drop of rate in such genes has a higher probability of fixation. α c YGR192CYOR248WYGR165WYGL223CYBR200WYEL077CYMR031CYNL095CYJR053WYER106W
Fig. 7. Reduction in steady-state output rate (production rate) R in percent as a function of α c ∈ [10 − , − ] for S. cerevisiae genes.
V. D
ISCUSSION
In many important processes biological “particles” move along some kind of a one-dimensional “track”.Examples include gene transcription and translation, cellular transport, and more. The flow can be eitherbidirectional (as in the case of transcription) or unidirectional (as in the case of translation), with thepossibility of both attachment and detachment of particles at different sites along the track. For example,motor proteins like kinesin and dynein that move along a certain microtubule may detach and attach toan overlapping microtubule.To rigorously model and analyze such processes, we introduced a new deterministic mathematical modelthat can be derived as the dynamic mean-field approximation of ASEP with Langmuir kinetics, called theMFALK. Our main results show that the MFALK is a monotone and contractive dynamical system. Thisimplies that it admits a globally asymptotically unique equilibrium point, and that it entrains to periodicexcitations (with a common period
T > ) in any of its rates, i.e. the densities along the chain, as wellas the output rate, converge to unique period solutions with period T .It is important to note that several known models are special cases of the MFALK. These include forexample the RFM [48], the model used in [15] for DNA transcription, and the model of phosphorelaysin [11]. Topics for further research include the following. In the RFM, it has been shown that the steady-state Although in this model the occupancy levels are normalized differently. production rate is related to the maximal eigenvalue of a certain non-negative, symmetric tridiagonalmatrix with elements that are functions of the RFM rates, i.e. the λ i s [44]. This implies that the mapping ( λ , . . . , λ n ) → R is strictly concave, and that sensitivity analysis of R is an eigenvalue sensitivity prob-lem [45]. An interesting research topic is whether R = R ( λ , . . . , λ n , γ , . . . , γ n , α , . . . , α n , β , . . . , β n ) in the MFALK can also be described using such a linear-algebraic approach.The application of the MFALK to model ribosome drop off suggests an interesting direction for furtherstudy, namely, how to design genes that minimize the drop off rate.Another research direction is motivated by the fact that many of the transport phenomena that canbe modeled using the MFALK do not take place in isolation. For example, many mRNA molecules aretranslated in parallel in the cell. Thus, a natural next step is to study networks of interconnected MFALKs.Graph theory can be used to describe the interconnections between the various MFALKs in the network.In this context, ribosome drop off may perhaps increase the total production rate in the entire system,as it allows ribosomes to detach from slow sites, enter the pool of free ribosomes, and then attach tothe initiation sites of other, less crowded, mRNA molecules. However, drop off still incurs the biological“cost” associated to the synthesis of a chain of amino-acids that is only a part of the desired protein. Thefact that the MFALK is contractive may prove useful in analyzing networks of MFALKs, as there existinteresting results proving the overall contractivity of a network based on contractivity of the subsystemsand their couplings (see, e.g. [5], [50]).Another interesting topic for further research is studying the effect of controlled detachment rates onthe formation of traffic jams. Indeed, it is known that kinesin-family motor proteins are more susceptibleto dissociation when their pathway is blocked [14], [62].A PPENDIX : P
ROOFS
We begin by discussing some symmetry properties of the MFALK, as these will be useful in the proofslater on. Symmetry
The MFALK enjoys two symmetries that will be useful later on. First, let z i ( t ) := 1 − x i ( t ) , i = 1 , . . . , n .In other words, z i ( t ) is the amount of “free space” at site i at time t . Then using (1) yields ˙ z = γ (1 − z ) + λ z (1 − z ) + α (1 − z ) − γ z (1 − z ) − λ z − β z , ˙ z = γ z (1 − z ) + λ z (1 − z ) + α (1 − z ) − γ z (1 − z ) − λ z (1 − z ) − β z , ... ˙ z n = γ n − z n − (1 − z n ) + λ n (1 − z n ) + α n (1 − z n ) − γ n z n − λ n − z n (1 − z n − ) − β n z n . (27)This is just the MFALK (1), but with the parameters permuted as follows: λ k → γ k , γ k → λ k , β k → α k ,and α k → β k for all k . The symmetry here follows from the fact that we can replace the roles of theforward and backward flows in the MFALK.Next, let y i ( t ) := 1 − x n +1 − i ( t ) , i = 1 , . . . , n . In other words, y i ( t ) is the amount of “free space” atsite n + 1 − i at time t . Then using (1) yields ˙ y = λ n (1 − y ) + γ n − y (1 − y ) + α n (1 − y ) − λ n − y (1 − y ) − γ n y − β n y , ˙ y = λ n − y (1 − y ) + γ n − y (1 − y ) + α n − (1 − y ) − λ n − y (1 − y ) − γ n − y (1 − y ) − β n − y , ... ˙ y n = λ y n − (1 − y n ) + γ (1 − y n ) + α (1 − y n ) − λ y n − γ y n (1 − y n − ) − β y n . (28)This is just the MFALK (1), but with the parameters permuted as follows: λ k → λ n − k , γ k → γ n − k , β k → α n +1 − k , and α k → β n +1 − k for all k . Note that (27) is simply (28) with the variable renaming z i → y n +1 − i , i = 1 , . . . , n .Both symmetries are reminiscent of the particle-hole symmetry in ASEP [6], [30]: the basic idea is thatthe progression of a particle from left to right is also the progression of a hole from right to left.ProofofProposition1. If (4) holds then the MFALK satisfies property ( BR ) in [35], and [35, Lemma 1]implies (6). If (5) holds then (27) satisfies property ( BR ) in [35], and this implies (6).Proof of Proposition 2. Write the MFALK as ˙ x = f ( x ) . A calculation shows that the Jacobian ma-trix J ( x ) := ∂f∂x ( x ) satisfies J ( x ) = L ( x ) + P , where L ( x ) is given in (30), and P is the diagonal L ( x ) = − λ (1 − x ) − γ x λ x + γ (1 − x ) . . . λ (1 − x ) + γ x − λ x − γ (1 − x ) − λ (1 − x ) − γ x . . . λ (1 − x ) + γ x . . . ... . . .
00 0 . . . λ n − x n − + γ n − (1 − x n − )0 0 . . . − λ n − x n − − γ n − (1 − x n − ) (30) matrix P = diag( − λ − γ − α − β , − α − β , . . . , − α n − − β n − , λ n − γ n − α n − β n ) . (29)Note that L ( x ) is tridiagonal and Metzler (i.e, every off-diagonal entry is non-negative) for any x ∈ C n .Recall that the matrix measure µ : R n × n → R induced by the L norm is given by µ ( A ) =max { c ( A ) , . . . , c n ( A ) } , where c i ( A ) is the sum of the elements in column i of A with off-diagonalelements taken with absolute value [64]. For the Jacobian J of the MFALK, µ ( J ( x )) = η for all x ∈ C n .It is well-known (see, e.g., [3]) that this implies (7).ProofofProposition3. For ζ ∈ [0 , / , let C nζ := { x ∈ C n : ζ ≤ x i ≤ − ζ , i = 1 , . . . , n } . Note that C n = C n , and that C nζ is a strict subcube of C n for all ζ ∈ (0 , / . By Proposition 1, forany τ > there exists ζ = ζ ( τ ) ∈ (0 , / , with ζ ( τ ) → as τ → , such that x ( t + τ, a ) ∈ C nζ , for all t ≥ and all a ∈ C n . (31)For any x ∈ C nζ every entry L ij on the sub- and super-diagonal of L in (30) satisfies L ij ≥ ζ s , where s :=min ≤ i ≤ n − { λ i + γ i } > . Combining this with [35, Theorem 4], implies that for any ζ ∈ (0 , / thereexists ε = ε ( ζ ) > , and a diagonal matrix D = diag(1 , q , q q , . . . , q q . . . q n − ) , with q i = q i ( ε ) > , such that the MFALK is contractive on C nζ w.r.t. the scaled L norm defined by | z | ,D := | Dz | .Furthermore, we can choose ε such that ε ( ζ ) → as ζ → , and D ( ε ) → I as ε → . Now Thm. 1in [39] implies that the MFALK is contractive after a small overshoot and short transient (SOST). Prop. 4in [39] implies that for the MFALK SOST is equivalent to SO, and this completes the proof. ProofofProposition4. We begin by recursively defining two sequences. For all integers i ≥ , let u i +1 = 1 + ℓ + ℓ + · · · + ℓ i ,ℓ i +1 = u i + ℓ + ℓ + · · · + ℓ i − . (32)with initial conditions u = u = 1 , and ℓ = 0 , ℓ = 1 . We claim that for k = 0 , , . . . , n − thesteady-state density in site n − k is generically the ratio of two polynomials in R : e n − k = p k ( R ) q k ( R ) , with deg( p k ( R )) = u k , deg( q k ( R )) = ℓ k . (33)We prove this by induction on k . By (16), e n = aR + b , with a := ( λ n + γ n + β n + α n ) − and b := ( γ n + β n ) a ,and this proves (33) for k = 0 . Using (16) again yields e n − = R + γ n − e n λ n − (1 − e n ) + γ n − e n = R + γ n − ( aR + b ) λ n − + ( γ n − − λ n − )( aR + b ) , and this proves (33) for k = 1 . Now assume that there exists s ≥ such that (33) holds for k =0 , , . . . , s − . By (16), e n − s = R + γ n − s e n − s +1 − g n − s +1 ( e n − s +1 ) − g n − s +2 ( e n − s +2 ) − · · · − g n − ( e n − ) λ n − s (1 − e n − s +1 ) + γ n − s e n − s +1 , and applying (12) and the induction hypothesis yields e n − s = R + γ n − s p s − q s − + ( β n − s +1 + α n − s +1 ) p s − q s − + ( β n − s +2 + α n − s +2 ) p s − q s − + · · · + ( β n − + α n − ) p q + cλ n − s + ( γ n − s − λ n − s ) p s − q s − , where c := − β n − s +1 −· · ·− β n − . Multiplying the numerator and the denominator by q . . . q s − yields e n − s = p s /q s , where deg( p s ) = max { q . . . q s − ) , deg( p s − q . . . q s − ) , . . . , deg( p q . . . q s − ) } , deg( q s ) = max { deg( q . . . q s − ) , deg( p s − q . . . q s − ) } . By the induction hypothesis, deg( p s ) = max { ℓ + · · · + ℓ s − , u s − + ℓ + · · · + ℓ s − , . . . , u + ℓ + · · · + ℓ s − } , deg( q s ) = max { ℓ + · · · + ℓ s − , u s − + ℓ + · · · + ℓ s − } . (34)It is straightforward to prove that (32) implies that ℓ i ≤ u i ≤ ℓ i + 1 , i = 0 , , , . . . . (35)Combining this with (34) yields deg( p s ) = 1 + ℓ + · · · + ℓ s − , and deg( q s ) = u s − + ℓ + · · · + ℓ s − .Thus, deg( p s ) = u s and deg( q s ) = ℓ s , and this completes the inductive proof of (33). In particular, (33)yields e = p n − ( R ) q n − ( R ) , (36)with deg( p n − ( R )) = u n − , deg( q n − ( R )) = ℓ n − . Substituting this in the last equation of (16) yields v p n − q n − = z − R + n − X j =2 g j ( e j ) , where v := λ + γ + β + α , and z := λ + β . Arguing as above shows that this is a polynomialequation of the form w ( R ) = 0 , with deg( w ) = 1 + ℓ + · · · + ℓ n − = u n . It is straightforward to proveby induction that (32) implies that u k = 1 + (cid:22) k (cid:23) , ℓ k = 2 k − ( − k , (we note in passing that the latter sequence is known as the Jacobsthal sequence [58]), and this completesthe proof of Proposition 4.Proof of Proposition 5. We begin by proving that R > implies that Q ni =0 λ i > Q ni =0 γ i . If R > then (18) yields λ (1 − e ) > γ e ,λ i e i (1 − e i +1 ) > γ i e i +1 (1 − e i ) , i = 1 , . . . , n − ,λ n e n > γ n (1 − e n ) . (37) Multiplying all these inequalities, and using the fact that e ∈ int( C n ) yields n Y i =0 λ i > n Y i =0 γ i . (38)To prove the converse implication, assume that (38) holds. Multiplying both sides of this inequality bythe strictly positive term Q nj =1 e j (1 − e j ) yields n Y i =0 a i > n Y i =0 b i , where a := λ (1 − e ) , a i := λ i e i (1 − e i +1 ) , i = 1 , . . . , n − , a n = λ n e n , b := γ e , b i := γ i e i +1 (1 − e i ) , i = 1 , . . . , n − , and b n = γ n (1 − e n ) . This means that a ℓ > b ℓ for some index ℓ ∈ { , . . . , n } . Since R = a ℓ − b ℓ (see (18)), it follows that R > . Summarizing, we showed that R > if and onlyif Q ni =0 λ i > Q ni =0 γ i . The proof that R < if and only if Q ni =0 λ i < Q ni =0 γ i is similar. This impliesthat R = 0 if and only if Q ni =0 λ i = Q ni =0 γ i . This completes the proof of (19).To prove (20), note that (18) yields e n = R + γ n λ n + γ n ,e i = R + γ i e i +1 λ i (1 − e i +1 ) + γ i e i +1 , i = n − , . . . , ,e = λ − Rλ + γ . (39)Substituting R = 0 completes the proof of Prop. 5.ProofofProposition6. Since the Jacobian J ( x ) of the MFALK is Metzler (i.e, every off-diagonal entryis non-negative) for any x ∈ C n , the MFALK is a cooperative system [60], and this yields (24).When λ i + γ i > , i = 1 , . . . , n − , the matrix L ( x ) and, therefore, J ( x ) , is irreducible for every x ∈ int( C n ) , and combining this with Proposition 1 implies (25) (see, e.g., [60, Ch. 4]).Proof of Theorem 1. The Jacobian of the PMFALK is J ( t, x ( t )) = L ( t, x ( t )) + P ( t ) , with L givenin (30), and P is given in (29) (but now with time-varying rates). Pick an initial time t ≥ , and τ > .The stated conditions guarantee the existence of ζ ∈ (0 , / such that x ( t, t , a ) ∈ C nζ for all t ≥ t + τ and all a ∈ C n . Also, [35, Thm. 4] implies that there exists a diagonally-scaled L norm such thatthe PMFALK is contractive on C nζ w.r.t. this norm. Now entrainment follows from known results oncontractive systems with a periodic excitation (see, e.g. [49]).ProofofProposition7. First, using Remark 1 and the argument used in the proof of [46, Prop. 4] shows that all the derivatives in the statement of of Proposition 7 exist.Given a MFALK, pick j ∈ { , . . . , n } and consider the new MFALK obtained by changing α j to ˜ α j ,with ˜ α j > α j , and all other rates unchanged. Let ˜ e , ˜ R denote the steady-state density and production ratein the modified MFALK. Seeking a contradiction, assume that ˜ e n ≥ e n . (40)Then (14) implies that ˜ R ≥ R, (41)and if j = n then ˜ R > R . By (15) with i = n − , R = λ n − e n − (1 − e n ) − γ n − e n (1 − e n − ) and ˜ R = λ n − ˜ e n − (1 − ˜ e n ) − γ n − ˜ e n (1 − ˜ e n − ) , and combining this with (40) and (41) yields ˜ e n − ≥ e n − . (42)Now using (15) with i = n − yields ˜ e n − ≥ e n − , and ˜ e n − > e n − if j = n − . Proceeding in this wayshows that ˜ e k ≥ e k , k = n, n − , . . . , j, (43) ˜ e k > e k , k = j − , j − , . . . , . (44)Combining this with (15) with i = 0 yields ˜ R < R . This contradicts (41), so ˜ e n > e n . (45)Proceeding as above yields ˜ e i > e i for all i , so ∂e i ∂α j < for all i, j . The proofs of all the other equationsin Prop. 7 are very similar and therefore omitted.R EFERENCES [1] B. Alberts, A. Johnson, J. Lewis, M. Raff, K. Roberts, and P. Walter,
Molecular Biology of the Cell . New York: Garland Science,2008.[2] R. J. R. Algar, T. Ellis, and G. Stan, “Modelling essential interactions between synthetic genes and their chassis cell,” in
Proc. 53rdIEEE Conference Decision and Control , LA, 2014, pp. 5437–5444.[3] Z. Aminzare and E. D. Sontag, “Contraction methods for nonlinear systems: A brief introduction and some open problems,” in
Proc.53rd IEEE Conf. on Decision and Control , Los Angeles, CA, 2014, pp. 3835–3847. [4] Y. Arava, Y. Wang, J. D. Storey, C. L. Liu, P. O. Brown, and D. Herschlag, “Genome-wide analysis of mRNA translation profiles inSaccharomyces cerevisiae,” Proceedings of the National Academy of Sciences , vol. 100, no. 7, pp. 3889–3894, 2003.[5] M. Arcak, “Certifying spatially uniform behavior in reaction-diffusion PDE and compartmental ODE systems,”
Automatica , vol. 47,no. 6, pp. 1219–1229, 2011.[6] R. A. Blythe and M. R. Evans, “Nonequilibrium steady states of matrix-product form: a solver’s guide,”
J. Phys. A: Math. Gen. , vol. 40,no. 46, pp. R333–R441, 2007.[7] Y. Chadani, K. Ono, S. Ozawa, Y. Takahashy, K. Takay, H. Nanamiya, Y. Tozawa, K. Kutsukake, and T. Abo, “Ribosome rescue byEscherichia coli arf-a (yhdl) in the absence of trans-translation systems,”
Mol. Microbiol. , vol. 78, pp. 796–808, 2010.[8] A. C. M. Cheung and P. Cramer, “Structural basis of RNA polymerase II backtracking, arrest and reactivation,”
Nature , vol. 471, no.7337, pp. 249–253, 03 2011.[9] D. Chu, E. Kazana, N. Bellanger, T. Singh, M. F. Tuite, and T. von der Haar, “Translation elongation can control translation initiationon eukaryotic mRNAs,”
EMBO J. , vol. 33, no. 1, pp. 21–34, 2014.[10] L. S. Churchman and J. S. Weissman, “Nascent transcript sequencing visualizes transcription at nucleotide resolution,”
Nature , vol.469, no. 7330, pp. 368–373, 01 2011.[11] A. Csikasz-Nagy, L. Cardelli, and O. S. Soyer, “Response dynamics of phosphorelays suggest their potential utility in cell signaling,”
J. Royal Society Interface , vol. 8, pp. 480–488, 2011.[12] A. Dana and T. Tuller, “Mean of the typical decoding rates: a new translation efficiency index based on the analysis of ribosomeprofiling data,” G3 , vol. 5, no. 1, pp. 73–80, 2014.[13] C. Desoer and H. Haneda, “The measure of a matrix as a tool to analyze computer algorithms for circuit analysis,” IEEE Trans. CircuitTheory , vol. 19, pp. 480–486, 1972.[14] R. Dixit, J. L. Ross, Y. E. Goldman, and E. L. Holzbaur, “Differential regulation of dynein and kinesin motor proteins by tau,”
Science ,vol. 319, pp. 1086–1089, 2008.[15] S. Edri, E. Gazit, E. Cohen, and T. Tuller, “The RNA polymerase flow model of gene transcription,”
IEEE Trans. Biomedical Circuitsand Systems , vol. 8, no. 1, pp. 54–64, 2014.[16] T. Ezaki and K. Nishinari, “Exact stationary distribution of an asymmetric simple exclusion process with Langmuir kinetics and memoryreservoirs,”
Journal of Physics A: Mathematical and Theoretical , vol. 45, no. 18, p. 185002, 2012.[17] M. J. Garca-Meseguer, J. A. V. de Labra, M. Garcia-Moreno, F. Garca-Canovas, B. H. Havsteen, and R. Varon, “Mean residence timesin linear compartmental systems. Symbolic formulae for their direct evaluation,”
Bull. Math. Biol. , vol. 65, pp. 279–308, 2003.[18] M. A. Gilchrist and A. Wagner, “A model of protein translation including codon bias, nonsense errors, and ribosome recycling,”
J.Theoretical Biology , vol. 239, no. 4, pp. 417–34, 2006.[19] F. S. Heldt, C. A. Brackley, C. Grebogi, and M. Thiel, “Community control in cellular protein production: consequences for amino acidstarvation,”
Philosophical Trans. Royal Society of London A: Mathematical, Physical and Engineering Sciences , vol. 373, no. 2056,2015.[20] S. D. Hooper and O. Berg, “Gradients in nucleotide and codon usage along Escherichia coli genes,”
Nucleic Acids Res , vol. 28, pp.3517–3523, 2000.[21] J. A. Jacquez and C. P. Simon, “Qualitative theory of compartmental systems,”
SIAM Review , vol. 35, no. 1, pp. 43–79, 1993.[22] F. Jorgensen and C. Kurland, “Processivity errors of gene expression in Escherichia coli,”
J. Mol. Biol. , vol. 215, pp. 511–521, 1990.[23] T. V. Karpinets, D. J. Greenwood, C. E. Sams, and J. T. Ammons, “RNA: protein ratio of the unicellular organism as a characteristic ofphosphorous and nitrogen stoichiometry and of the cellular requirement of ribosomes for protein synthesis,”
BMC Biol. , vol. 4, no. 30,pp. 274–80, 2006. [24] K. Keiler, “Mechanisms of ribosome rescue in bacteria,” Nat. Rev. Microbiol. , vol. 13, pp. 285–297, 2015.[25] K. Keiler, P. Waller, and R. Sauer, “Role of a peptide tagging system in degradation of proteins synthesised from damaged messengerrna,”
Science , vol. 271, pp. 990–993, 1996.[26] B. S. Kerner, “The physics of green-wave breakdown in a city,”
Europhysics Letters , vol. 102, no. 2, p. 28010, 2013.[27] S. Klumpp and T. Hwa, “Traffic patrol in the transcription of ribosomal RNA,”
RNA Biol. , vol. 6, no. 4, pp. 392–4, 2009.[28] C. Kurland, “Translational accuracy and the fitness of bacteria,”
Ann. Rev. Genet. , vol. 26, pp. 29–50, 1992.[29] C. Kurland and R. Mikkola,
Starvation in Bacteria . NY: Plenum Press, 1993, ch. The impact of nutritional state on the microevolutionof ribosomes, pp. 225–238.[30] G. Lakatos, T. Chou, and A. Kolomeisky, “Steady-state properties of a totally asymmetric exclusion process with periodic structure,”
Phys. Rev. E , vol. 71, p. 011103, 2005.[31] C. Leduc, K. Padberg-Gehle, V. Varga, D. Helbing, S. Diez, and J. Howard, “Molecular crowding creates traffic jams of kinesin motorson microtubules,”
Proceedings of the National Academy of Sciences , vol. 109, pp. 6100–6105, 2012.[32] R. Lipowsky and S. Klumpp, “Life is motion: multiscale motility of molecular motors,”
Physica A: Statistical Mechanics and itsApplications , vol. 352, no. 1, pp. 53–112, 2005.[33] R. Lipowsky, S. Klumpp, and T. M. Nieuwenhuizen, “Random walks of cytoskeletal motors in open and closed compartments,”
Phys.Rev. Lett. , vol. 87, no. 10, p. 108101, 2001.[34] W. Lohmiller and J.-J. E. Slotine, “On contraction analysis for non-linear systems,”
Automatica , vol. 34, pp. 683–696, 1998.[35] M. Margaliot, E. D. Sontag, and T. Tuller, “Entrainment to periodic initiation and transition rates in a computational model for genetranslation,”
PLoS ONE , vol. 9, no. 5, p. e96039, 2014.[36] M. Margaliot and T. Tuller, “On the steady-state distribution in the homogeneous ribosome flow model,”
IEEE/ACM Trans.Computational Biology and Bioinformatics , vol. 9, pp. 1724–1736, 2012.[37] M. Margaliot and T. Tuller, “Stability analysis of the ribosome flow model,”
IEEE/ACM Trans. Computational Biology andBioinformatics , vol. 9, pp. 1545–1552, 2012.[38] M. Margaliot and T. Tuller, “Ribosome flow model with positive feedback,”
J. Royal Society Interface , vol. 10, p. 20130267, 2013.[39] M. Margaliot, E. D. Sontag, and T. Tuller, “Contraction after small transients,”
Automatica , vol. 67, pp. 178–184, 2016.[40] J. Mierczynski, “A class of strongly cooperative systems without compactness,”
Colloq. Math. , vol. 62, pp. 43–47, 1991.[41] E. Nudler, “RNA polymerase backtracking in gene regulation and genome instability,”
Cell , vol. 149, no. 7, pp. 1438–1445, 2015.[42] A. Parmeggiani, T. Franosch, and E. Frey, “Phase coexistence in driven one-dimensional transport,”
Phys. Rev. Lett. , vol. 90, p. 086601,2003.[43] A. Parmeggiani, T. Franosch, and E. Frey, “Totally asymmetric simple exclusion process with Langmuir kinetics,”
Phys. Rev. E , vol. 70,p. 046101, 2004.[44] G. Poker, Y. Zarai, M. Margaliot, and T. Tuller, “Maximizing protein translation rate in the nonhomogeneous ribosome flow model: aconvex optimization approach,”
J. Royal Society Interface , vol. 11, no. 100, 2014.[45] G. Poker, M. Margaliot, and T. Tuller, “Sensitivity of mRNA translation,”
Sci. Rep. , vol. 5, p. 12795, 2015.[46] A. Raveh, M. Margaliot, E. D. Sontag, and T. Tuller, “A model for competition for ribosomes in the cell,”
J. Royal Society Interface ,vol. 116, no. 20151062, 2016.[47] A. Raveh, Y. Zarai, M. Margaliot, and T. Tuller, “Ribosome flow model on a ring,”
IEEE/ACM Trans. Computational Biology andBioinformatics , vol. 12, no. 6, pp. 1429–1439, 2015.[48] S. Reuveni, I. Meilijson, M. Kupiec, E. Ruppin, and T. Tuller, “Genome-scale analysis of translation elongation with a ribosome flowmodel,”
PLOS Computational Biology , vol. 7, p. e1002127, 2011. [49] G. Russo, M. di Bernardo, and E. D. Sontag, “Global entrainment of transcriptional systems to periodic inputs,” PLOS ComputationalBiology , vol. 6, p. e1000739, 2010.[50] G. Russo, M. di Bernardo, and E. D. Sontag, “A contraction approach to the hierarchical analysis and design of networked systems,”
IEEE Trans. Automat. Control , vol. 58, no. 5, pp. 1328–1331, 2013.[51] I. W. Sandberg, “On the mathematical foundations of compartmental analysis in biology, medicine, and ecology,”
IEEE Trans. Circuitsand Systems , vol. 25, no. 5, pp. 273–279, 1978.[52] A. Schadschneider, D. Chowdhury, and K. Nishinari,
Stochastic Transport in Complex Systems: From Molecules to Vehicles . Elsevier,2011.[53] M. Schliwa and G. Woehlke, “Molecular motors,”
Nature , vol. 422, pp. 759–765, 2003.[54] J. W. Shaevitz, E. A. Abbondanzieri, R. Landick, and S. M. Block, “Backtracking by single RNA polymerase molecules observed atnear-base-pair resolution,”
Nature , vol. 426, no. 6967, pp. 684–687, 12 2003.[55] L. B. Shaw, R. K. P. Zia, and K. H. Lee, “Totally asymmetric exclusion process with extended objects: a model for protein synthesis,”
Phys. Rev. E , vol. 68, p. 021910, 2003.[56] C. Shoemaker, D. Eyler, and R. Green, “Dom34:hbs1 promotes subunit dissociation and peptidyl-tRNA drop-off to initiate no-godecay,”
Science , vol. 330, no. 6002, pp. 369–72, 2010.[57] C. Sin, D. Chiarugi, and A. Valleriani, “Quantitative assessment of ribosome drop-off in E. coli,”
Nucleic Acids Res. , vol. 44, no. 6,pp. 2528–37, 2016.[58] N. J. A. Sloane and S. Plouffe,
The Encyclopedia of Integer Sequences . Academic Press, 1995.[59] J. Smillie, “Competitive and cooperative tridiagonal systems of differential equations,”
SIAM J. Mathematical Analysis , vol. 15, pp.530–534, 1984.[60] H. L. Smith,
Monotone Dynamical Systems: An Introduction to the Theory of Competitive and Cooperative Systems , ser. MathematicalSurveys and Monographs. Providence, RI: Amer. Math. Soc., 1995, vol. 41.[61] A. Subramaniam, B. Zid, and E. O’Shea, “An integrated approach reveals regulatory controls on bacterial translation elongation,”
Cell ,vol. 159, no. 5, pp. 1200–11, 2014.[62] I. A. Telley, P. Bieling, and T. Surrey, “Obstacles on the microtubule reduce the processivity of kinesin-1 in a minimal in vitro systemand in cell extract,”
Biophys. J. , vol. 96, pp. 3341–3353, 2009.[63] T. Tuller and H. Zur, “Multiple roles of the coding sequence 5’ end in gene expression regulation,”
Nucleic Acids Res. , vol. 43, no. 1,pp. 13–28, 2015.[64] M. Vidyasagar,
Nonlinear Systems Analysis . Englewood Cliffs, NJ: Prentice Hall, 1978.[65] R. Willmann, G. Sch¨utz, and D. Challet, “Exact Hurst exponent and crossover behavior in a limit order market model,”
Physica A:Statistical Mechanics and its Applications , vol. 316, no. 1, pp. 430–440, 2002.[66] A. Yonath, “Ribosomes: Ribozymes that survived evolution pressures but is paralyzed by tiny antibiotics,” in
MacromolecularCrystallography: Deciphering the Structure, Function and Dynamics of Biological Molecules , A. M. Carrondo and P. Spadon, Eds.Dordrecht: Springer Netherlands, 2012, pp. 195–208.[67] Z. Zafrir, H. Zur, and T. Tuller, “Selection for reduced translation costs at the intronic 5’ end in fungi,”
DNA Res. , vol. 23, no. 4, pp.377–94, 2016.[68] H. Zaher and R. Green, “A primary role for elastase factor 3 in quality control during translation elongation in Escherichia coli,”
Cell ,vol. 147, pp. 396–408, 2011.[69] S. Zaher and R. Green, “Quality control by the ribosome following peptide bond formation,”
Nature , vol. 457, pp. 161–166, 2009.[70] Y. Zarai, M. Margaliot, E. D. Sontag, and T. Tuller, “Controlling mRNA translation,” 2016, submitted. [71] Y. Zarai, M. Margaliot, and T. Tuller, “Explicit expression for the steady-state translation rate in the infinite-dimensional homogeneousribosome flow model,” IEEE/ACM Trans. Computational Biology and Bioinformatics , vol. 10, pp. 1322–1328, 2013.[72] G. Zhang, I. Fedyunin, O. Miekley, A. Valleriani, A. Moura, and Z. Ignatova, “Global and local depletion of ternary complex limitstranslational elongation,”
Nucleic Acids Res. , vol. 38, no. 14, pp. 4778–4787, 2010.[73] R. K. P. Zia, J. Dong, and B. Schmittmann, “Modeling translation in protein synthesis with TASEP: A tutorial and recent developments,”
J. Statistical Physics , vol. 144, pp. 405–428, 2011.[74] A. Zupanic, C. Meplan, S. M. Grellscheid, J. C. Mathers, T. B. Kirkwood, J. E. Hesketh, and D. P. Shanley, “Detecting translationalregulation by change point analysis of ribosome profiling data sets,”