A Model for Competition for Ribosomes in the Cell
Alon Raveh, Michael Margaliot, Eduardo D. Sontag, Tamir Tuller
aa r X i v : . [ q - b i o . GN ] A ug A Model for Competition for Ribosomes in the Cell
Alon Raveh, Michael Margaliot, Eduardo D. Sontag* and Tamir Tuller*
Abstract
A single mammalian cell includes an order of − mRNA molecules and as many as − ribosomes.Large-scale simultaneous mRNA translation and the resulting competition for the available ribosomes has importantimplications to the cell’s functioning and evolution. Developing a better understanding of the intricate correlationsbetween these simultaneous processes, rather than focusing on the translation of a single isolated transcript, shouldhelp in gaining a better understanding of mRNA translation regulation and the way elongation rates affect organismalfitness. A model of simultaneous translation is specifically important when dealing with highly expressed genes,as these consume more resources. In addition, such a model can lead to more accurate predictions that are neededin the interconnection of translational modules in synthetic biology.We develop and analyze a general model for large-scale simultaneous mRNA translation and competition forribosomes. This is based on combining several ribosome flow models (RFMs) interconnected via a pool of freeribosomes. We prove that the compound system always converges to a steady-state and that it always entrainsor phase locks to periodically time-varying transition rates in any of the mRNA molecules. We use this modelto explore the interactions between the various mRNA molecules and ribosomes at steady-state. We show thatincreasing the length of an mRNA molecule decreases the production rate of all the mRNAs. Increasing any of thecodon translation rates in a specific mRNA molecule yields a local effect: an increase in the translation rate of thismRNA, and also a global effect: the translation rates in the other mRNA molecules all increase or all decrease.These results suggest that the effect of codon decoding rates of endogenous and heterologous mRNAs on proteinproduction is more complicated than previously thought. Index Terms mRNA translation, competition for resources, systems biology, monotone dynamical systems, first integral,entrainment, synthetic biology, context-dependence in mRNA translation, heterologous gene expression.
I. I
NTRODUCTION
Various processes in the cell utilize the same finite pool of available resources. This means that theprocesses actually compete for these resources, leading to an indirect coupling between the processes.This is particularly relevant when many identical intracellular processes, all using the same resources,take place in parallel.Biological evidence suggests that the competition for RNA polymerase (RNAP) and ribosomes, andvarious transcription and translation factors, is a key factor in the cellular economy of gene expression.The limited availability of these resources is one of the reasons why the levels of genes, mRNA, andproteins produced in the cell do not necessarily correlate [63], [32], [57], [67], [53], [69], [29].It was estimated that in a yeast cell there are approximately , mRNA molecules. These can betranslated in parallel [74], [71], with possibly many ribosomes scanning the same transcript concurrently. The research of MM and TT is partially supported by a research grant from the Israeli Ministry of Science, Technology, and Space. Thework of EDS is supported in part by grants AFOSR FA9550-14-1-0060 and ONR 5710003367.*Corresponding authors: EDS and TT.A. Raveh 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]. D. Sontag is with the Dept. of Mathematics and the Center for Quantitative Biology, Rutgers University, Piscataway, NJ 08854, USA.E-mail: [email protected]. Tuller is with the Dept. of Biomedical Engineering and the Sagol School of Neuroscience, Tel-Aviv University, Tel-Aviv 69978, Israel.E-mail: [email protected]
Fig. 1. Illustration of simultaneous translation of mRNA chains (right) interconnected via a pool of free ribosomes (left).
The number of ribosomes is limited (in a yeast cell it is approximately , ) and this leads to acompetition for ribosomes. For example, if more ribosomes bind to a certain mRNA molecule then thepool of free ribosomes in the cell is depleted, and this may lead to lower initiation rates in the othermRNAs (see Fig. 1).There is a growing interest in computational or mathematical models that take into account the com-petition for available resources in translation and/or transcription (see, for example, [22], [23], [44], [9],[14], [7]). One such model, that explicitly considers the movement of the ribosomes [RNAP] along themRNA [DNA], is based on a set of asymmetric simple exclusion processes (ASEPs) interconnected to apool of free ribosomes. ASEP is an important model from non-equilibrium statistical physics describingparticles that hop randomly from one site to the next along an ordered lattice of sites, but only if thenext site is empty. This form of “rough exclusion” models the fact that the particles cannot overtake oneanother. ASEP has been used to model and analyze numerous multiagent systems with local interactionsincluding the flow of ribosomes along the mRNA molecule [38],[58]. In this context, the lattice representsthe mRNA molecule, and the particles are the ribosomes. For more on mathematical and computationalmodels of translation, see the survey paper [70].Ref. [24] considered a closed system composed of a single ASEP connected to a pool (or reservoir) of“free” particles. The total number of particles is conserved. This is sometimes referred to as the parkinggarage problem , with the lattice modeling a road, the particles are cars, and the pool corresponds to aparking garage. Ref. [12] studied a similar system using domain wall theory. Ref. [13] (see also [21])considered a network composed of two ASEPs connected to a finite pool of particles. The analysis inthese papers focuses on the phase diagram of the compound system with respect to certain parameters,and on how the phase of one ASEP affects the phase of the other ASEPs. These studies rely on the phasediagram of a single ASEP that is well-understood only in the case where all the transition rates inside thechain (the elongation rates) are equal. Thus, the network is typically composed of homogeneous ASEPs.Another model [7] combines non-homogenous ASEPs in order to study competition between multiplespecies of mRNA molecules for a pool of tRNA molecules. This study was based on the
Saccharomyces cerevisiae genome. However, in this case (and similar models, such as [14]) analysis seems intractableand one must resort to simulations only.Our approach is based on the ribosome flow model (RFM) [52]. This is a deterministic, continuous-time,synchronous model for translation that can be derived via the mean-field approximation of ASEP [6].The RFM includes n state-variables describing the ribosomal density in n consecutive sites along themRNA molecule, and n + 1 positive parameters: the initiation rate λ , and the elongation rate λ i fromsite i to site i + 1 , i = 1 , . . . , n . The RFM has a unique equilibrium point e = e ( λ , . . . , λ n ) , and anytrajectory emanating from a feasible initial condition converges to e [42] (see also [39]). This means thatthe system always converges to a steady-state ribosomal density that depends on the rates, but not on theinitial condition. In particular, the production rate converges to a steady-state value R . The mapping fromthe rates to R is a concave function, so maximizing R subject to a suitable constraint on the rates is aconvex optimization problem [48], [73]. Sensitivity analysis of the RFM with respect to the rates has beenstudied in [49]. These results are important in the context of optimizing the protein production rate insynthetic biology. Ref. [39] has shown that when the rates λ i are time-periodic functions, with a commonminimal period T , then every state-variable converges to a periodic solution with period T . In otherwords, the ribosomal densities entrain to periodic excitations in the rates (due e.g. to periodically-varyingabundances of tRNA molecules).In ASEP with periodic boundary conditions a particle that hops from the last site returns to the first one.The mean field approximation of this model is called the ribosome flow model on a ring (RFMR). Theperiodic boundary conditions mean that the total number of ribosomes is conserved. Ref. [51] analyzedthe RFMR using the theory of monotone dynamical systems that admit a first integral.Both the RFM and the RFMR model mRNA translation on a single mRNA molecule. In this paper,we introduce a new model, called the RFM network with a pool (RFMNP), that includes a networkof RFMs, interconnected through a dynamical pool of free ribosomes, to model and analyze simultaneoustranslation and competition for ribosomes in the cell. To the best of our knowledge, this is the first studyof a network of RFMs. The total number of ribosomes in the RFMNP is conserved, leading to a firstintegral of the dynamics. Applying the theory of monotone dynamical systems that admit a first integralwe prove several mathematical properties of the RFMNP: it admits a continuum of equilibrium points,every trajectory converges to an equilibrium point, and any two solutions emanating from initial conditionscorresponding to an equal number of ribosomes converge to the same equilibrium point. These resultshold for any set of rates and in particular when the RFMs in the network are not necessarily homogeneous.These stability results are important because they provide a rigorous framework for studying questionssuch as how does a change in one RFM affects the behavior of all the other RFMs in the network? Indeed,since a steady-state exists, this can be reduced to asking how does a change in one RFM in the networkaffects the steady-state behavior of the network?To analyze competition for ribosomes, we consider the effect of increasing one of the rates in one RFM,say RFM λ i in RFM tends to increase the steady-statedensity in sites i + 1 , i + 2 , . . . , and decrease the density in site i of this RFM. The total density (i.e., thesum of all the densities on the different sites along RFM ) can either decrease or increase. In the firstcase, more ribosomes are freed to the pool, and this increases the initiation rates in all the other RFMsleading to higher production rates. The second case leads to the opposite result. The exact outcome ofincreasing one of the rates thus depends on the many parameter values defining the pool and the set ofRFMs in the network.Our model takes into account the dynamics of the translation elongation stage, yet is still amenableto analysis. This allows to develop a rigorous understanding of the effect of competition for ribosomes.Previous studies on this topic were either based on simulations (see, for example, [14], [7]) or did notinclude a dynamical model of translation elongation (see, e.g., [23], [9]). For example, in an interesting paper, combining mathematical modeling and biological experiments, Gyorgy et al. [23] study the expres-sion levels of two adjacent reporter genes on a plasmid in E. coli based on measurements of fluorescencelevels, that is, protein levels. These are of course the result of all the gene expression steps (transcription,translation, mRNA degradation, protein degradation) making it difficult to separately study the effect ofcompetition for ribosomes or to study specifically the translation elongation step. Their analysis yieldsthat the attainable output p , p of the two proteins satisfies the formula αp + βp = Y, (1)where Y is related to the total number of ribosomes (but also other translation factors and possiblyadditional gene expression factors), and α, β are constants that depend on parameters such as the plasmidcopy number, dissociation constants of the ribosomes binding to the Ribosomal Binding Site (RBS), etc.This equation implies that increasing the production of one protein always leads to a decrease in theproduction of the other protein (although more subtle correlations may take place via the effects on theconstants α and β ). A similar conclusion has been derived for other models as well [16, Ch. 7].In our model, improving the translation rate of a codon in one mRNA may either increase or decreasethe translation rates of all other mRNAs in the cell. Indeed, the effect on the other genes depends on thechange in the total density of ribosomes on the modified mRNA molecule, highlighting the importanceof modeling the dynamics of the translation elongation step. We show however that when increasing theinitiation rate in an RFM in the network, the total density in this RFM always increases and, consequently,the production rate in all the other RFMs decreases. This special case agrees with the results in [23].Another recent study [50] showed that the hidden layer of interactions among genes arising fromcompetition for shared resources can dramatically change network behavior. For example, a cascade ofactivators can behave like an effective repressor, and a repression cascade can become bistable. This agreeswith several previous studies in the field (see, for example, [29], [2]).The remainder of this paper is organized as follows. The next section describes the new model. Wedemonstrate using several examples how it can be used to study translation at the cell level. Section IIIdescribes our main theoretical results, and details their biological implications. To streamline the presen-tation, the proofs of the results are placed in the Appendix. The final section concludes and describesseveral directions for further research.II. T HE MODEL AND SOME EXAMPLES
Since our model is based on a network of interconnected RFMs, we begin with a brief review ofthe RFM.
A. Ribosome flow model
The RFM models the traffic flow of ribosomes along the mRNA. The mRNA chain is divided into aset of n compartments or sites, where each site may correspond to a codon or a group of codons. Thestate-variable x i ( t ) , i = 1 , . . . , n , describes the ribosome occupancy at site i at time t , normalized suchthat x i ( t ) = 0 [ x i ( t ) = 1 ] implies that site i is completely empty [full] at time t . Roughly speaking, onemay also view x i ( t ) as the probability that site i is occupied at time t . The dynamical equations of theRFM are: ˙ x = λ (1 − x ) − λ x (1 − x ) , ˙ x = λ x (1 − x ) − λ x (1 − x ) , ˙ x = λ x (1 − x ) − λ x (1 − x ) , ... ˙ x n − = λ n − x n − (1 − x n − ) − λ n − x n − (1 − x n ) , ˙ x n = λ n − x n − (1 − x n ) − λ n x n . (2) R ( t ) = λ n x n ( t ) ProductionProtein λ λ λ λ λ n − λ n x ( t ) x ( t ) x ( t ) x n ( t ) CodonSite
Fig. 2. Topology of the RFM. State variable x i ( t ) ∈ [0 , describes the normalized ribosome occupancy level in site i at time t . Theinitiation rate is λ , and λ i is the elongation rate between sites i and i + 1 . Production rate at time t is R ( t ) := λ n x n ( t ) . These equations describe the movement of ribosomes along the mRNA chain. The transition rates λ , . . . , λ n are all positive numbers (units=1/time). To explain this model, consider the equation ˙ x = λ x (1 − x ) − λ x (1 − x ) . The term λ x (1 − x ) represents the flow of particles from site to site . This isproportional to the occupancy x at site and also to − x , i.e. the flow decreases as site becomesfuller. In particular, if x = 1 , i.e. site is completely full, the flow from site to site is zero. This isa “soft” version of the rough exclusion principle in ASEP. Note that the maximal possible flow rate fromsite to site is the transition rate λ . The term λ x (1 − x ) represents the flow of particles from site to site .The dynamical equations for the other state-variables are similar. Note that λ controls the initiationrate into the chain, and that R ( t ) := λ n x n ( t ) , is the rate of flow of ribosomes out of the chain, that is, the translation (or protein production) rate attime t . The RFM topology is depicted in Fig. 2.The RFM encapsulates simple exclusion and unidirectional movement along the lattice just as in ASEP.This is not surprising, as the RFM can be derived via a mean-field approximation of ASEP (see, e.g., [6,p. R345] and [62, p. 1919]). However, the analysis of these two models is quite different as the RFMis a deterministic, continuous-time, synchronous model, whereas ASEP is a stochastic, discrete-type,asynchronous one.In order to study a network of interconnected RFMs, it is useful to first extend the RFM into a single-input single-output (SISO) control system: ˙ x = λ (1 − x ) u − λ x (1 − x ) , ˙ x = λ x (1 − x ) − λ x (1 − x ) , ˙ x = λ x (1 − x ) − λ x (1 − x ) , ... ˙ x n − = λ n − x n − (1 − x n − ) − λ n − x n − (1 − x n ) , ˙ x n = λ n − x n − (1 − x n ) − λ n x n ,y = λ n x n . (3)Here the translation rate becomes the output y of the system, and the flow into site is multiplied bya time-varying control u : R + → R + , representing the flow of ribosomes from the “outside world” intothe strand which is related to the rate ribosomes diffuse to the 5’end (in eukaryotes) or the RBS (inprokaryotes) of the mRNA. Of course, mathematically one can absorb λ into u , but we do not do thisbecause we think of λ as representing some local/mRNA-specific features of the mRNA sequence (e.g.the strength of the Kozak sequences in eukaryotes or the RBS in prokaryote).The set of admissible controls U is the set of bounded and measurable functions taking values in R + for all time t ≥ . Eq. (3), referred to as the RFM with input and output (RFMIO) [43], facilitates thestudy of RFMs with feedback connections. We note in passing that (3) is a monotone control system asdefined in [4]. From now on we write (3) as ˙ x = f ( x, u ) ,y = λ n x n . (4)Let C n := { z ∈ R n : z i ∈ [0 , , i = 1 , . . . , n } , denote the closed unit cube in R n . Since the state-variables represent normalized occupancy levels, wealways consider initial conditions x (0) ∈ C n . It is straightforward to verify that C n is an invariant setof (4), i.e. for any u ∈ U and any x (0) ∈ C n the trajectory satisfies x ( t, u ) ∈ C n for all t ≥ . B. RFM network with a pool
To model competition for ribosomes in the cell, we consider a set of m ≥ RFMIOs, represent-ing m different mRNA chains. The i th RFMIO has length n i , input function u i , output function y i , andrates λ i , . . . , λ in i . The RFMIOs are interconnected through a pool of free ribosomes (i.e., ribosomes thatare not attached to any mRNA molecule). The output of each RFMIO is fed into the pool, and the poolfeeds the initiation locations in the mRNAs (see Fig. 3). Thus, the model includes m RFMIOs: ˙ x = f ( x , u ) , y = λ n x n , ... ˙ x m = f ( x m , u m ) , y m = λ mn m x n m , (5)and a dynamic pool of ribosomes described by ˙ z = m X j =1 y j − m X j =1 λ j (1 − x j ) G j ( z ) . (6)where z : R + → R describes the pool occupancy. Eq. (6) means that the flow into the pool is the sumof all output rates P mj =1 y j of the RFMIOs minus the total flow of ribosomes that bind to an mRNAmolecule P mj =1 λ j (1 − x j ) G j ( z ) . Recall that the term (1 − x j ) represents the exclusion, i.e. as the firstsite in RFMIO j becomes fuller, less ribosomes can bind to it. Thus, the input to RFMIO j is u j ( t ) = G j ( z ( t )) , j = 1 , . . . , m. (7)We assume that each G j ( · ) : R + → R + satisfies: (1) G j (0) = 0 ; G j is continuously differentiableand G ′ j ( z ) > for all z ≥ (so G j is strictly increasing on R + ); and (3) there exists s > suchthat G j ( z ) ≤ sz for all z > sufficiently small. These properties imply in particular that if the pool isempty then no ribosomes can bind to the mRNA chains, and that as the pool becomes fuller the initiationrates to the RFMIOs increase.Typical examples for functions satisfying these properties include the linear function, say, G j ( z ) = z ,and G j ( z ) = a j tanh( b j z ) , with a j , b j > . In the first case, the flow of ribosomes into the first siteof RFM i is given by λ i z (1 − x i ) , and the product here can be justified via mass-action kinetics. Theuse of tanh may be suitable for modeling a saturating function. This is in fact a standard function in ASEPmodels with a pool [13], [1], because it is zero when z is zero, uniformly bounded, and strictly increasingfor z ≥ . Also, for z ≥ the function tanh( z ) takes values in [0 , so it can also be interpreted as aprobability function [20].In the context of a shared pool, it is natural to consider the special case where G j ( z ) = G ( z ) forall j = 1 , . . . , m . The differences between the initiation sites in the strands are then modeled by thedifferent λ j ’s. RFMIO
P OOL y m y y u u m u G m ( z ) G ( z ) G ( z ) Fig. 3. Topology of the RFMNP.
Note that combining the properties of G j with (6) implies that if z (0) ≥ then z ( t ) ≥ for all t ≥ .Thus, the pool occupancy is always non-negative.Summarizing, the RFM network with a pool (RFMNP) is given by equations (5), (6), and (7). This isa dynamical system with d := 1 + P mi =1 n i state-variables. Example 1
Consider a network with m = 2 RFMIOs, the first [second] with dimension n = 2 [ n = 3 ].Then the RFMNP is given by ˙ x = λ (1 − x ) G ( z ) − λ x (1 − x ) , ˙ x = λ x (1 − x ) − λ x , ˙ x = λ (1 − x ) G ( z ) − λ x (1 − x ) , ˙ x = λ x (1 − x ) − λ x (1 − x ) , ˙ x = λ x (1 − x ) − λ x , ˙ z = λ x + λ x − λ (1 − x ) G ( z ) − λ (1 − x ) G ( z ) . (8)Note that this system has d = 6 state-variables. (cid:3) An important property of the RFMNP is, that being a closed system, the total occupancy H ( t ) := z ( t ) + m X j =1 n j X i =1 x ji ( t ) , (9)is conserved, that is, H ( t ) ≡ H (0) , for all t ≥ . (10)In other words, H is a first integral of the dynamics. In particular, this means that z ( t ) ≤ H ( t ) = H (0) for all t ≥ , i.e. the pool occupancy is uniformly bounded.The RFMNP models mRNAs that compete for ribosomes because the total number of ribosomes isconserved. As more ribosomes bind to the RFMs, the pool depletes, G j ( z ) decreases, and the effectiveinitiation rate to all the RFMs decreases. This allows to systematically address important biologicalquestions on large-scale simultaneous translation under competition for ribosomes. The following examples demonstrate this. We prove in Section III that all the state-variables in the RFMNP converge to a steady-state. Let e ij ∈ [0 , denote the steady-state occupancy in site j in RFMIO i , and let e z ∈ [0 , ∞ ) denotethe steady-state occupancy in the pool. In the examples below we always consider these steady-statevalues (obtained numerically by simulating the differential equations). Example 2
Although we are mainly interested in modeling large-scale simultaneous translation, it isnatural to first consider a model with a single mRNA molecule connected to a pool of ribosomes. From abiological perspective, this models the case where there is one gene that is highly expressed with respectto all other genes (e.g. an extremely highly expressed heterologous gene).Consider an RFMNP that includes a single RFMIO (i.e. m = 1 ), with dimension n = 3 , rates λ i = 1 , i = 0 , , , , and a pool with output function G ( z ) = tanh( z ) . We simulated this system for the initialcondition x i (0) = 0 , z (0) = c for various values of c . Note that H ( t ) ≡ H (0) = c . Fig. 4 depicts thesteady-state values e , e , e of the state-variables in the RFMIO, and the steady-state pool occupancy e z .It may be seen that for small values of c the steady-state ribosomal densities and thus the productionrates are very low. This is simply because there are not enough ribosomes in the network. The ribosomaldensities increase with c . For large values of c , the output function of the pool saturates, as tanh( z ) → ,and so does the initiation rate in the RFMIO. Thus, the densities in the RFMIO saturate to the valuescorresponding to the initiation rate λ = 1 , and then all the remaining ribosomes accumulate in the pool.Using a different pool output function, for example G ( z ) = z , leads to the same qualitative behavior, butwith higher saturation values for the ribosomal densities in the RFMIO. (Note that the ribosomal densitiesin an RFM are finite even when λ → ∞ [41].) (cid:3) This simple example already demonstrates the coupling between the ribosomal pool, initiation rate,and elongation rates. When the ribosomal pool is small the initiation rate is low. Thus, the ribosomaldensities on the mRNA are low and there are no interactions between ribosomes (i.e., no “traffic jams”)along the mRNA. The initiation rate becomes the rate limiting step of translation. On the other-hand,when there are many ribosomes in the pool the initiation rate increases, the elongation rates become ratelimiting and “traffic jams” along the mRNA evolve. At some point, a further increase in the number ofribosomes in the pool will have a negligible effect on the production rate.It is known that there can be very large changes in the number of ribosomes in the cell during e.g.exponential growth. For example. Ref. [8] reports changes in the range , to , . The exampleabove demonstrates how these large changes in the number of ribosomes are expected to affect thetranslational regimes; specifically, it may cause a switch between the different regimes mentioned above.The next example describes an RFMNP with several mRNA chains. Let n ∈ R n denote the vectorof n ones. Example 3
Consider an RFMNP with m = 3 RFMIOs of dimensions n = n = n = 3 , and rates λ i = c, λ i = 5 , λ i = 10 , i = 0 , . . . , . In other words, every RFMIO has homogeneous rates. Suppose also that G i ( z ) = tanh( z ) , for i = 1 , , .We simulated this RFMNP for different values of c with the initial condition z (0) = 0 , x (0) = (1 / , x (0) = (1 / and x (0) = (1 / . Thus, H (0) = 3 . in all the simulations. For each value of c ,every state-variable in the RFMNP converges to a steady-state. Fig. 5 depicts the steady-state value e z and the steady-state output y i in each RFMIO. It may be seen that increasing c , i.e. increasing all theelongation rates in RFMIO leads to an increase in the steady-state translation rates in all the RFMIOsin the network. Also, it leads to an increase in the steady-state occupancy of the pool. It may seem thatthis contradicts (10) but this is not so. Increasing c indeed increases all the steady-state translation rates,but it decreases the steady-state occupancies inside each RFMIO so that the total H ( t ) = H (0) = 3 . isconserved.Define the average steady-state occupancy (ASSO) in RFMIO j by ¯ e j := n j P n j i =1 e ji . Fig. 6 depictsthe ASSO in each RFMIO as a function of c . It may be seen as c increases the ASSO in RFMIO c e e e e z Fig. 4. Steady-state values in the RFMNP in Example 2 as a function of the total occupancy H (0) = c . c y y y z Fig. 5. Steady-state outputs y i of the three RFMIOs and pool occupancy z as a function of the homogeneous transition rate c in RFMIO . decreases quickly, yet the ASSOs in the other two RFMIOs slowly increases. Since the ribosomes spendless time on RFMIO (due to increased c ) they are now available for translating the other RFMIOs,leading to the increased ASSO in the other mRNAs. (cid:3) From a biological point of view this example corresponds to a situation where accelerating one ofthe mRNA chains increases the protein production rates in all the mRNAs and also increases the numberof free ribosomes. Surprisingly, perhaps, it also suggests that a relatively larger number of free ribosomesin the cell corresponds to higher protein production rates. This agrees with evolutionary, biological, andsynthetic biology studies that have suggested that (specifically) highly expressed genes (that are transcribedinto many mRNA molecules) undergo selection to include codons with improved elongation rates [32],[66], [63]. Specifically, two mechanisms by which improved codons affect translation efficiency and theorganismal fitness are [66]: (1) global mechanism : selection for improved codons contributes towardimproved ribosomal recycling and global allocation; the increased number of free ribosomes improvesthe effective translation initiation rate of all genes, and thus improves global translation efficiency; and(2) local mechanism : the improved translation elongation rate of an mRNA contributes directly to itsprotein production rate. c e e e Fig. 6. Average steady-state occupancy in the three RFMIOs as a function of the homogeneous transition rate c in RFMIO . The example above demonstrates both mechanisms, as improvement of the translation elongation ratesof one RFM increases the translation rate of this mRNA (local translation efficiency), and also of the otherRFMs (global translation efficiency). In addition, as can be seen, the decrease in ASSO in RFMIO issignificantly higher than the increase in ASSO in the other RFMIO. Thus, the simulation also demonstratesthat increasing the translation rate c may contribute to decreasing ribosomal collision (and possiblyribosomal abortion).We prove in Section III that when one of the rates in one of the RFMIOs increases two outcomesare possible: either all the production rates in the other RFMs increase (as in this example) or they alldecrease. As discussed below, we believe that this second case is less likely to occur in endogenous genes,but may occur in heterologous gene expression.The next example describes the effect of changing the length of one RFMIO in the network. Example 4
Consider an RFMNP with m = 2 RFMIOs of dimensions n and n = 10 , rates λ i = 1 , i = 0 , . . . , n ,λ j = 1 , j = 0 , . . . , , and G i ( z ) = tanh( z/ , i = 1 , . In other words, both RFMIOs have the same homogeneous rates.We simulated this RFMNP for different values of n with the initial condition z (0) = 100 , x (0) = 0 n , x (0) = 0 . Thus H (0) = 100 in all the simulations. For each value of n , every state-variable inthe RFMNP converges to a steady-state. Fig. 7 depicts the steady-state values of z , and the steady-stateoutput y i in each RFMIO. It may be seen that increasing n , i.e. increasing the length of RFMIO leads to a decrease in the steady-state production rates and in the steady-state pool occupancy. This isreasonable, as increasing n means that ribosomes that bind to the first chain remain on it for a longerperiod of time. This decreases the production rate y and, by the competition for ribosomes, also decreasesthe pool occupancy and thus decreases y . (cid:3) From a biological point of view this suggests that decreasing the length of mRNA molecules contributeslocally and globally to improving translation efficiency. A shorter coding sequence improves the translationrate of the mRNA and, by competition, may also improve the translation rates in all other mRNAs. Thus,we should expect to see selection for shorter coding sequences, specifically in highly expressed genes andin organisms with large population size. Indeed, previous studies have reported that in some organismsthe coding regions of highly expressed genes tend to be shorter [14]; other studies have shown that other n y y z/300 Fig. 7. Steady-state outputs y i of the two RFMIOs and steady-state pool occupancy z/ as a function of the length n of RFMIO .The normalization of z by is used to obtain similar magnitudes for all the values only. (non-coding) parts of highly expressed genes tend to be shorter [11], [17], [35], [64]. Decreasing thelength of different parts of the gene should contribute to organismal fitness via improving the energeticcost of various gene expression steps. For example, shorter genes should improve the metabolic cost ofsynthesizing mRNA and proteins; it can also reduce the energy spent for splicing and processing of RNAand proteins. However, there are of course various functional and regulatory constraints that also contributeto shaping the gene length (see, for example [10]). Our results and these previous studies suggest thatin some cases genes are expected to undergo selection also for short coding regions, as this reduces therequired number of translating ribosomes.The next section describes theoretical results on the RFMNP. All the proofs are placed in the Appendix.III. M ATHEMATICAL PROPERTIES OF THE
RFMNPLet
Ω := [0 , n × · · · × [0 , n m × [0 , ∞ ) denote the state-space of the RFMNP (recall that every x ji takes values in [0 , and z ∈ [0 , ∞ ) ). Foran initial condition a ∈ Ω , let (cid:2) x ( t, a ) z ( t, a ) (cid:3) ′ denote the solution of the RFMNP at time t . It isstraightforward to show that the solution remains in Ω for all t ≥ . Our first result shows that a slightlystronger property holds. A. Persistence
Proposition 1
For any τ > there exists ε = ε ( τ ) > , with ε ( τ ) → when τ → , such that forall t ≥ τ , all j = 1 , . . . , m , all i = 1 , . . . , n j , and all a ∈ (Ω \ { } ) , ε ≤ x ji ( t, a ) ≤ − ε, and ε ≤ z ( t, a ) . In other words, after time τ the solution is ε -separated from the boundary of Ω . This result is usefulbecause on the boundary of Ω , denoted ∂ Ω , the RFMNP looses some desirable properties. For example,its Jacobian matrix may become reducible on ∂ Ω . Prop. 1 allows us to overcome this technical difficulty,as it implies that any trajectory is separated from the boundary after an arbitrarily short time. B. Strong Monotonicity
Recall that a cone K ⊆ R n defines a partial order in R n as follows. For two vectors a, b ∈ R n , wewrite a ≤ b if ( b − a ) ∈ K ; a < b if a ≤ b and a = b ; and a ≪ b if ( b − a ) ∈ Int( K ) . A dynamicalsystem ˙ x = f ( x ) is called monotone if a ≤ b implies that x ( t, a ) ≤ x ( t, b ) for all t ≥ . In other words,monotonicity means that the flow preserves the partial ordering [59]. It is called strongly monotone if a < b implies that x ( t, a ) ≪ x ( 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 order is called cooperative .The next result analyzes the cooperativity of the RFMNP. Let d := 1 + m P i =1 n i denote the dimension ofthe RFMNP. Proposition 2
For any a, b ∈ Ω with a ≤ b , x ( t, a ) ≤ x ( t, b ) and z ( t, a ) ≤ z ( t, b ) , for all t ≥ . (11) Furthermore, if a < b then x ( t, a ) ≪ x ( t, b ) and z ( t, a ) < z ( t, b ) , for all t > . (12)This means the following. Consider the RFMNP initiated with two initial conditions such that theribosomal densities in every site and the pool corresponding to the first initial condition are smaller orequal to the densities in the second initial condition. Then this correspondence between the densitiesremains true for all time t ≥ . C. Stability
For s ≥ , let L s := { y ∈ Ω : 1 ′ d y = s } . In other words, L s is a level set of the first integral H . Theorem 1
Every level set L s , s ≥ , contains a unique equilibrium point e L s of the RFMNP, and for anyinitial condition a ∈ L s , the solution of the RFMNP converges to e L s . Furthermore, for any ≤ s < p , e L s ≪ e L p . (13)In particular, this means that every trajectory converges to an an equilibrium point, representing steady-state ribosomal densities in the RFMIOs and the pool. Note that Proposition 1 implies that for any s > , e L s ∈ Int(Ω) . Eq. (13) means that the continuum of equilibrium points, namely, { e L s : s ∈ [0 , ∞ ) } ,are linearly ordered. Example 5
Consider an RFMNP with m = 2 RFMIOs with dimensions n = n = 1 , and G i ( z ) = z , i = 1 , , i.e. ˙ x = λ (1 − x ) z − λ x , ˙ x = λ (1 − x ) z − λ x , ˙ z = λ x + λ x − λ (1 − x ) z − λ (1 − x ) z. (14)Note that even in this simple case the RFMNP is a nonlinear system. Assume that λ = λ = 1 , andthat λ = λ , and denote this value simply by λ . Pick an initial condition in Ω , and let s := x (0)+ x (0)+ x x z Fig. 8. Trajectories of the system in Example 5 for three initial conditions in L . The equilibrium point (15) is marked by a circle. z (0) , so that the trajectory belong to L s for all t ≥ . Any equilibrium point e = (cid:2) e e e z (cid:3) ′ ∈ L s satisfies (1 − e ) e z = λe , (1 − e ) e z = λe ,e + e + e z = s. This yields two solutions e = e = (cid:16) s + 2 + λ − p ( s + 2 + λ ) − s (cid:17) / , (15) e z = (cid:16) s − − λ + p ( s + 2 + λ ) − s (cid:17) / , and e = e = (cid:16) s + 2 + λ + p ( s + 2 + λ ) − s (cid:17) / ,e z = (cid:16) s − − λ − p ( s + 2 + λ ) − s (cid:17) / , It is straightforward to verify that in the latter solution e z < , so this is not a feasible solution. Thesolution (15) does belong to L s , so the system admits a unique equilibrium in L s . Fig. 8 depicts trajectoriesof (14) for three initial conditions in L , namely, (cid:2) (cid:3) , (cid:2) (cid:3) , and (cid:2) / / (cid:3) , and theequilibrium point (15) for s = λ = 1 . It may be seen that every one of these trajectories converges to e . (cid:3) D. Contraction
Contraction theory is a powerful tool for analyzing nonlinear dynamical systems (see, e.g., [36]), withapplications to many models from systems biology [3], [54], [40]. In a contractive system, the distancebetween any two trajectories decreases at an exponential rate. It is clear that the RFMNP is not a contractivesystem on Ω , with respect to any norm, as it admits more than a single equilibrium point. Nevertheless,the next result shows that the RFMNP is non-expanding with respect to the ℓ norm: | q | = P di =1 | q i | . Proposition 3
For any a, b ∈ Ω , (cid:12)(cid:12)(cid:12)(cid:12)(cid:20) x ( t, a ) z ( t, a ) (cid:21) − (cid:20) x ( t, b ) z ( t, b ) (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) ≤ | a − b | , for all t ≥ . (16)In other words, the ℓ distance between trajectories can never increase.Pick a ∈ Ω , and let s := 1 ′ d a . Substituting b = e L s in (16) yields (cid:12)(cid:12)(cid:12)(cid:12)(cid:20) x ( t, a ) z ( t, a ) (cid:21) − e L s (cid:12)(cid:12)(cid:12)(cid:12) ≤ | a − e L s | , for all t ≥ . (17)This means that the convergence to the equilibrium point e L s is monotone in the sense that the ℓ distance to e L s can never increase. Combining (17) with Theorem 1 implies that every equilibrium pointof the RFMNP is semistable [26]. E. Entrainment
Many important biological processes are periodic. Examples include circadian clocks and the cell cycledivision process. Proper functioning requires certain biological systems to follow these periodic patterns,i.e. to entrain to the periodic excitation.In the context of translation, it has been shown that both the RFM [39] and the RFMR [51] entrain toperiodic translation rates, i.e. if all the transition rates are periodic time-varying functions, with a common(minimal) period
T > then each state variable converges to a periodic trajectory, with a period T . Herewe show that the same property holds for the RFMNP.We say that a function f is T -periodic if f ( t + T ) = f ( t ) for all t . Assume that the λ ji ’s in the RFMNPare time-varying functions satisfying: • there exist < δ < δ such that λ ji ( t ) ∈ [ δ , δ ] for all t ≥ and all j ∈ { , . . . , m } , i ∈ { , . . . , n j } . • there exists a (minimal) T > such that all the λ ji ’s are T -periodic.We refer to the model in this case as the periodic ribosome flow model network with a pool (PRFMNP). Theorem 2
Consider the PRFMNP. Fix an arbitrary s > . There exists a unique function φ s : R + → Int(Ω) , that is T -periodic, and for any a ∈ L s the solution of the PRFMNP converges to φ s . In other words, every level set L s of H contains a unique periodic solution, and every solution ofthe PRFMNP emanating from L s converges to this solution. Thus, the PRFMNP entrains (or phaselocks) to the periodic excitation in the λ ji ’s. This implies in particular that all the protein production ratesconverge to a periodic pattern with period T .Note that since a constant function is a periodic function for any T , Theorem 2 implies entrainment toa periodic trajectory in the particular case where one of the λ ji ’s oscillates, and all the other are constant.Note also that the stability result in Theorem 1 follows from Theorem 2. Example 6
Consider the RFMNP (8) with G i ( z ) = tanh( z ) , and all rates equal to one except for λ ( t ) =5 + 4 sin(2 πt ) . In other words, there is a single time-varying periodic rate in RFMIO . Note that allthese rates are periodic with a common minimal period T = 1 . Fig. 9 depicts the solution of this PRFMNPas a function of time t for . ≤ t ≤ . The initial condition is z (0) = x ij (0) = 1 / for all i, j . It maybe seen that all the state variables converge to a periodic solution. In particular, all state variables x i ( t ) converge to a periodic solution with (minimal) period T = 1 , and so does the pool occupancy z ( t ) .The x j ( t ) ’s also converge to a periodic solution, but it is not possible to tell from the figure whether thereare small oscillations with period T = 1 or the convergence is to a constant (of course, in both cases thisis a periodic solution with period T = 1 ). However, it can be shown using the first two equations in (8)that if z ( t ) converges to a periodic solution then so do x ( t ) and x ( t ) . Note that the peaks in x ( t ) arecorrelated with dips in x ( t ) , this is because when λ ( t ) is high on some time interval, i.e. the transition t x x x x x z Fig. 9. Trajectories of the PRFMNP in Example 6 as a function of time. rate from site to site is high, there is a high flow of ribosomes from site to site during thisinterval. (cid:3) From the biophysical point of view this means that the coupling between the mRNA molecules caninduce periodic oscillations in all the protein production rates even when all the transition rates in themolecules are constant, except for a single rate in a single molecule that oscillates periodically. Thetranslation rate of codons is affected among others by the tRNA supply (i.e. the intracellular abundanceof the different tRNA species) and demand (i.e. total number of codons from each type on all the mRNAmolecules) (see for example [47]). Thus, the translation rate of a codon(s) is affected by changes in thedemand (e.g. oscillations in mRNA levels) or by changes in the supply (e.g. oscillations in tRNA levels).The results reported here may suggest that oscillations in the mRNA levels of some genes or in theconcentration of some tRNA species (that occur for example during the cell cycle [60], [19]), can induceoscillations in the translation rates of the rest of the genes.
F. Competition
We already know that any trajectory of the RFMNP converges to an equilibrium point. A naturalquestion is how will a change in the parameters (that is, the transition rates) affect this equilibrium point.For example, if we increase some transition rate λ ji in RFMIO j , how will this affect the steady-stateproduction rate in the other RFMIOs? Without loss of generality, we assume that the change is in atransition rate of RFMIO . Theorem 3
Consider an RFMNP with m RFMIOs with dimensions n , . . . , n m . Let λ := (cid:2) λ . . . λ mn m (cid:3) ′ denote the set of all parameters of the RFMNP, and let e = (cid:2) e . . . e n e . . . e n . . . e m . . . e mn m e z (cid:3) ′ ∈ (0 , n + ··· + n m × R ++ denote the equilibrium point of the RFMNP on some fixed level set of H . Pick i ∈ { , . . . , n } . Considerthe RFMNP obtained by modifying λ i to ¯ λ i , with ¯ λ i > λ i . Let ¯ e denote the equilibrium point in the newRFMNP and let ˜ e := ¯ e − e . Then ˜ e i < , and ˜ e j > , for all j ∈ { i + 1 , ..., n } , (18) sign(˜ e ij ) = sign(˜ e z ) , for all i = 1 and all j. (19)(In the case i = 0 , the condition ˜ e i < is vacuous.) Increasing λ i means that ribosomes flow “more easily” from site i to site i + 1 in RFMIO . Eq. (18)means that the effect on the density in this RFMIO is that the number of ribosomes in site i decreases,and the number of ribosomes in all the sites to the right of site i increases. Eq. (19) describes the effecton the steady-state densities in all the other RFMIOs and the pool: either all these steady-state valuesincrease or they all decrease. The first case agrees with the results in Example 3 above.Note that (18) does not provide information on the change in e j , j < i . Our simulations show that anyof these values may either increase or decrease, with the outcome depending on the various parametervalues. Thus, the amount of information provided by (18) depends on i . In particular, when λ n is changedto ¯ λ n > λ n then the information provided by (18) is only that ˜ e n < . Much more information is available when i = 0 . Corollary 1
Suppose that λ is changed to ¯ λ > λ . Then ˜ e j > , for all j ∈ { , ..., n } , (20) and ˜ e ij < , for all i = 1 and all j, and ˜ e z < . (21)Indeed, for i = 0 , (18) yields (20). Also, we know that the changes in the densities in all other RFMIOsand the pool have the same sign. This sign cannot be positive, as combining this with (20) contradictsthe conservation of ribosomes, so (21) follows.In other words, increasing λ yields an increase in all the densities in RFM , and a decrease in allthe other densities. This makes sense, as increasing λ means that it is easier for ribosomes to bind to themRNA molecule. This increases the total number of ribosomes along this molecule and, by competition,decreases all the densities in the other molecules and the pool. Note that this special case agrees wellwith the results described in [23] (see (1)).It is important to emphasize, however, that there are various possible intracellular mechanisms that mayaffect λ ji , i > . For example, synonymous mutation/changes (in endogenous or heterologous) genes insidethe coding region may affect the adaptation of codons to the tRNA pool (codons that are recognized bytRNA with higher intracellular abundance usually tend to be translated more quickly [15]), the local foldingof the mRNA (stronger folding tend to decrease elongation rate [65]), or the interaction/hybridizationbetween the ribosomal RNA and the mRNA [34] (there are nucleotides sub-sequence that tend to interactwith the ribosomal RNA, causing transient pausing of the ribosome, and delay the translation elongationrate). Non synonymous mutation/changes inside the coding region may also affect the elongation forexample via the interaction between the nascent peptide and the exit tunnel of the ribosome [37], [55].In addition, intracellular changes in various translation factors (e.g. tRNA levels, translation elongationfactors, concentrations of amino acids, concentrations of Aminoacyl tRNA synthetase) and, as explainedabove, the mRNA levels can also affect elongation rates. Furthermore, various recent studies have demon-strated that manipulating the codons of a heterologous gene tend to result in significant changes in thetranslation rates and protein levels of the gene [32], [72], [5].Thus, our study is relevant to fundamental biological phenomena that are not covered in models thatdo not take into account the elongation dynamics. Example 7
Consider the RFMNP in (8) with G i ( z ) = z , λ = λ = 1 , λ = λ = 0 . , λ = 1 , and initialcondition (1 / . We consider a range of values for λ . For each fixed value, we simulated the dynamicsuntil steady-state for two cases: λ = 1 and ¯ λ = 10 . Fig. 10 depicts ˜ e , ˜ e for the various fixed valuesof λ . It may be seen that we always have ˜ e < and ˜ e > . Fig. 11 depicts ˜ e i , i = 1 , , , and ˜ e z forthe various fixed values of λ . It may be seen that for a small value of λ all the ˜ e i ’s and ˜ e z are negative,whereas for large values of λ they all become positive. (cid:3) λ ˜ e ˜ e Fig. 10. Steady-state perturbations in RFM in the RFMNP in Example 7 when changing λ = 1 to ¯ λ = 10 for various values of λ . −3 λ ˜ e ˜ e ˜ e ˜ e z Fig. 11. Steady-state perturbations in RFM in the RFMNP in Example 7 when changing λ = 1 to ¯ λ = 10 for various values of λ . Theorem 3 implies that when the codons of a gene are modified into “faster codons” (either via syntheticengineering or during evolution) then either all the translation rates of the other genes increase, or they alldecrease. However, Theorem 3 does not provide information on when each of these two cases happens.In order to address this, we need to calculate derivatives of the equilibrium point coordinates with respectto the rates. The next result shows that these derivatives are well-defined. Denote the mapping from theparameters to the unique equilibrium point in
Int(Ω) by α , that is, e ij = α ij ( λ, H (0)) , i = 1 , . . . , m , j = 1 , . . . , n i . Proposition 4
The derivative ∂∂λ pq α ij ( λ, H (0)) exists for all i, j, p, q . The next example uses these derivatives to obtain information on the two cases that can take place aswe change one of the rates.
Example 8
Consider an RFMNP with m = 2 RFMIO’s with lengths n and ℓ . To simplify the notation,let e = (cid:2) e , . . . , e n (cid:3) ′ [ v = (cid:2) v , . . . , e ℓ (cid:3) ′ ] denote the equilibrium point of RFMIO [RFMIO ], andlet λ i , i = 0 , . . . , n , denote the rates along RFMIO . Suppose that λ is changed to ¯ λ . Differentiating the steady-state equations λ G ( e z )(1 − e ) = λ e (1 − e ) = ... = λ n − e n − (1 − e n ) = λ n e n , n X i =1 e i + ℓ X j =1 v j + e z = H (0) , w.r.t. λ yields λ G ′ ( e z ) e ′ z (1 − e ) − λ G ( e z ) e ′ = λ n e ′ n ,e ′ z + ℓ X j =1 v ′ j = − n X i =1 e ′ i , where we use the notation f ′ := ∂∂λ f . These two equations yield ( λ n + λ G ′ ( e z )(1 − e )) e ′ z + λ n ℓ X j =1 v ′ j = ( λ G ( e z ) − λ n ) e ′ − λ n n − X i =2 e ′ i . Recall that G ( e z ) > , G ′ ( e z ) > , λ j > for all j , and < e p < for all p . Also, by Theorem 3, e ′ < , e ′ j > , for all j ∈ , ..., n , and sign( e ′ z ) = sign( v ′ k ) , for all k . Thus, sign( e ′ z ) = sign( v ′ j ) = sign ( λ G ( e z ) − λ n ) e ′ − λ n n − X i =2 e ′ i ! . (22)This means that the sign of the change in the densities in all the other RFMIO’s and the pool dependson several steady-state quantities including terms related to the initiation rate λ G ( e z ) and exit rate λ n in RFMIO , and also the change in the total density P n − i =1 e ′ i in this RFMIO.In the particular case n = 2 (i.e., a very short RFM), Eq. (22) becomes: sign( e ′ z ) = sign( v ′ j ) = sign( λ − λ G ( e z )) . (23)Note that λ G ( e z ) [ λ ] is the steady-state initiation [exit] rate in RFMIO . Thus, λ − λ G ( e z ) > means that it is “easier” for ribosomes to exit than to enter RFMIO , and in this case (23) means thatwhen λ is increased the change in all other densities will be positive. This is intuitive, as more ribosomeswill exit the modified molecule and this will improve the production rates in the other molecules. On theother-hand, if λ − λ G ( e z ) < then it is “easier” for ribosomes to enter than to exit RFMIO , soincreasing λ will lead to an increased number of ribosomes in RFMIO and, by competition, to adecrease in the production rate in all the other RFMIOs. (cid:3) Note that in the example above increasing λ always increases the steady-state production rate R = λ e in RFM (recall that e ′ > ). One may expect that this will always lead to an increase in the productionrate in the second RFMIO as well. However, the behavior in the RFMNP is more complicated becausethe shared pool generates a feedback connection between the RFMIO’s in the network. In particular, theeffect on the other RFMIO’s depends not only on the modified production rate of RFMIO , but onother factors including the change in the total ribosome density in RFMIO (see (22)).This analysis of a very short RFM suggests that the steady-state initiation rate of the mRNA with themodified codon plays an important role in determining the effect of modifications in the network. If thisinitiation rate is relatively low (so it becomes the rate limiting factor), as believed to be the case in mostendogenous genes [28], then the increase in the rate of one codon of the mRNA increases the translationrate in all the other mRNAs, whereas when this initiation rate is high then the opposite effect is obtained.This latter case may occur for example when a heterologous gene is highly expressed and thus “consumes”some of the available elongation/termination factors making the elongation rates the rate limiting factors. IV. D
ISCUSSION
We introduced a new model, the RFMNP, for large-scale simultaneous translation and competition forribosomes that combines several RFMIOs interconnected via a dynamic pool of free ribosomes. To the bestof our knowledge, this is the first model of a network composed of interconnected RFMIOs. The RFMNPis amenable to analysis because it is a monotone dynamical system that admits a non-trivial first integral.Our results show that the RFMNP has several nice properties: it is an irreducible cooperative dynamicalsystem admitting a continuum of linearly ordered equilibrium points, and every trajectory converges to anequilibrium point. The RFMNP is also on the “verge of contraction” with respect to the ℓ norm, and itentrains to periodic transition rates with a common period. The fact that the total number of ribosomes inthe network is conserved means that local properties of any mRNA molecule (e.g., the abundance of thecorresponding tRNA molecules) affects its own translation rate, and via competition, also globally affectsthe translation rates of all the other mRNAs in the network.An important implication of our analysis and simulation results is that there are regimes and parametervalues where there is a strong coupling between the different “translation components” (ribosomes andmRNAs) in the cell. Such regimes cannot be studied using models for translation of a single isolatedmRNA molecule. The RFMNP is specifically important when studying highly expressed genes with manymRNA molecules and ribosomes translating them because the dynamics of such genes strongly affectsthe ribosomal pool. For example, changes in the translation dynamics of a heterologous gene which isexpressed with a very strong promoter, resulting in very high mRNA copy number should affect theentire tRNA pool, and thus the translation of other endogenous genes. Highly expressed endogenousgenes ”consume” many ribosomes. Thus, a mutation that affects their (“local”) translation rate is expectedto affect also the translation dynamics of other mRNA molecules. Studying the evolution of such genesshould be based on understanding the global effect of such mutations using a computational model suchas the RFMNP.On the other hand, we can approximate the dynamics of genes that are not highly expressed (e.g., agene with mRNA levels that are . of the mRNA levels in the cell) using a single RFM. In this case,the relative effect of the mRNA on all other mRNAs is expected to be limited.Our analysis shows that increasing the translation initiation rate of a heterologous gene will alwayshave a negative effect on the translation rate of other genes (i.e. their translation rates decrease) and viceversa. The effect of increasing [decreasing] the translation rate of a codon of the heterologous gene on thetranslation rate of other genes is more complicated: while it always increases [decreases] the translationrate of the heterologous gene it may either increase or decrease the translation rate of all other genes. Thespecific outcome of such a manipulation can be studied using the RFMNP with the parameters based onthe heterologous genes and the host genome.Our analysis suggests that the effect of improving the transition rate of a codon in an mRNA moleculeon the production rate of other genes and the pool of ribosomes depends on the initiation rate in themodified mRNA. When the initiation rate is very low the effect is expected to be positive (all otherproduction rates increase). However, if the initiation rate is high the effect may be negative. This maypartially explain the selection for slower codons in highly expressed genes that practically decrease theinitiation rate [67], [63]. This relation may also suggest a new factor that contributes to the evolutionof highly expressed genes towards higher elongation and termination rates (i.e., the tendency of highlyexpressed genes to include “fast” codons). Indeed, lower elongation rates (and thus a relatively highinitiation rate) may decrease the production rates of other mRNAs that are needed for proper functioningof the organism.The RFM, and thus also the model described here, do not capture certain aspects of mRNA translation.For example, eukaryotic ribosomes may translate mRNAs in multiple cycles before entering the free ribo-somal pool [70], [43], [31]. This phenomenon may perhaps be modeled by adding positive feedback [43]in the RFMNP. In addition, different genes are transcribed at different rates, resulting in a different numberof (identical) mRNA copies for different genes. This can be modeled using a set of identical RFMs for each gene. Such a model can help in understanding how changes in mRNA levels of one gene affect thetranslation rates of all the mRNAs. The analysis here suggests that modifying the mRNA levels of a genewill affect the translation rates of all other genes in the same way. These and other aspects of biologicaltranslation may be integrated in our model in future studies. Ref. [25] develops the notion of the realizableregion for steady-state gene expression under resource limitations, and methods for mitigating the effectsof ribosome competition. Another interesting research direction is studying these topics in the context ofthe RFMNP.The results reported here can be studied experimentally by designing and expressing a library ofheterologous genes [72], [5], [32]. The effect of the manipulation of a codon (i.e, increasing or decreasingits rate) of the heterologous gene on the ribosomal densities and translation rates of all the mRNAs(endogenous and heterologous) can be performed via ribosome profiling [27] in addition to measurementsof mRNA levels, translation rates, and protein levels [56].We believe that networks of interconnected RFMIOs may prove to be powerful modeling and analysistools for other natural and artificial systems as well. These include communication networks, intracellulartrafficking in the cell, coordination of large groups of organisms (e.g., ants), traffic control, and more.A CKNOWLEDGMENTS
We thank Yoram Zarai for helpful comments.A
PPENDIX : P
ROOFS
ProofofProposition1. We require the following result on repelling boundaries and persistence.
Lemma 1
Consider a time-varying system ˙ x = f ( t, x ) (24) whose trajectories evolve on Ω := I × I × . . . × I n ⊆ R n + , where each I j is an interval of the form [0 , A ] , A > , or [0 , ∞ ) . Suppose that the time-dependent vector field f = (cid:2) f , . . . , f n (cid:3) ′ has the followingtwo properties: • the cyclic boundary-repelling property (CBR) : For each δ > and each sufficiently small ∆ > ,there exists K = K ( δ, ∆) > such that, for each k = 1 , . . . , n and each t ≥ , the condition x k ( t ) ≤ ∆ , and x k − ( t ) ≥ δ (25) (all indexes are modulo n ) implies that f k ( t, x ) ≥ K ; (26) • for any i ∈ { , . . . , n } , and any s ≥ , x i ( s ) > implies that x i ( t ) > , for all t ≥ s. (27) Then given any τ > there exists ε = ε ( τ ) > , with ε ( τ ) → as τ → , such that every solution x ( t, a ) ,with a = 0 , satisfies x i ( t, a ) ≥ ε for all i ∈ { , . . . , n } and all t ≥ τ. In other words, the conclusion is that after an arbitrarily short time every x i ( t, a ) is separated away fromzero.ProofofLemma1. Pick τ > and a = 0 . Then there exists i ∈ { , . . . , n } such that a i > . Since the( CBR ) condition is cyclic, we may assume without loss of generality that x n (0) > . Then (27) impliesthat there exists δ > such that x n ( t ) ≥ δ for all t ∈ [0 , τ /n ] .Fix ∆ > such that (CBR) holds. Let K = K ( δ, ∆) , and define ε := min { ∆ , Kτ /n } . Let t ∈ [0 , τ /n ] be such that x ( t ) ≥ ε . Such a t exists, since by property (CBR) , x ( t ) < ε ≤ ∆ for all t ∈ [0 , τ /n ] would imply that ˙ x ( t ) = f ( t, x ( t )) ≥ K > for all t ∈ [0 , τ /n ] , which in turn implies x ( τ /n ) ≥ x (0) + Kτ /n ≥ Kτ /n ≥ ε , contradicting x ( τ /n ) < ε . We claim that also x ( t ) ≥ ε forevery t ≥ t . Indeed, suppose otherwise. Then, there is some t > t such that ξ := x ( t ) < ε . Let σ := min { t ≥ t | x ( t ) ≤ ξ } > t . As x ( σ ) ≤ ξ < ε ≤ ∆ , and x n ( σ ) > property (CBR) says that ˙ x ( σ ) = f ( σ, x ( σ )) > , so it followsthat ˙ x ( t ) > on an interval [ σ − c, σ ] , for some c > . But then x ( σ − c ) < x ( σ ) ≤ ξ , which contradictsthe minimality of σ . Thus x ( t ) ≥ ε for all t ≥ t , and in particular x ( t ) ≥ ε , for all t ≥ τ /n. (28)Let ¯ K := K ( ε , ∆) , and define ε := min { ∆ , ¯ Kτ /n } . Let t ∈ [ τ /n, τ /n ] be such that x ( t ) ≥ ε .Such a t exists, since by property (CBR) and (28), x ( t ) < ε ≤ ∆ for all t ∈ [ τ /n, τ /n ] wouldimply that ˙ x ( t ) = f ( t, x ( t )) ≥ ¯ K > for all t ∈ [ τ /n, τ /n ] , which in turn implies x (2 τ /n ) ≥ x ( τ /n ) + ¯ Kτ /n ≥ ¯ Kτ /n ≥ ε , contradicting x (2 τ /n ) < ε . We claim that also x ( t ) ≥ ε for every t ≥ t . Indeed, suppose otherwise. Then, there is some t > t such that ξ := x ( t ) < ε . Let σ := min { t ≥ t | x ( t ) ≤ ξ } > t . As x ( σ ) ≤ ξ < ε ≤ ∆ , and x ( σ ) ≥ ε property (CBR) says that ˙ x ( σ ) = f ( σ, x ( σ )) > , so itfollows that ˙ x ( t ) > on an interval [ σ − c, σ ] , for some c > . But then x ( σ − c ) < x ( σ ) ≤ ξ , whichcontradicts the minimality of σ . Thus x ( t ) ≥ ε for all t ≥ t , and in particular x ( t ) ≥ ε , for all t ≥ τ /n. Continuing in this manner we have that for every i ∈ { , . . . , n } there exists ε i > such that x i ( t ) ≥ ε i , for all t ≥ iτ /n. Thus, x i ( t ) ≥ min { ε , . . . , ε n } , for all t ≥ τ and all i = 1 , . . . , n, and this completes the proof of Lemma 1. (cid:3) We can now prove Proposition 1. Consider the case m = 1 , i.e. the RFMNP is given by ˙ x = λ (1 − x ) G ( x n +1 ) − λ x (1 − x ) , ˙ x = λ x (1 − x ) − λ x (1 − x ) , ... ˙ x n − = λ n − x n − (1 − x n − ) − λ n − x n − (1 − x n ) , ˙ x n = λ n − x n − (1 − x n ) − λ n x n , ˙ x n +1 = λ n x n − λ (1 − x ) G ( x n +1 ) , (29)where we write x n +1 instead of z . The proof in the case where m > is similar. We begin by showingthat (29) satisfies the properties in Lemma 1 on Ω = [0 , n × [0 , ∞ ) . Fix an arbitrary δ > . If x ≤ ∆ and x n +1 ≥ δ then f = λ (1 − x ) G ( x n +1 ) − λ x (1 − x ) ≥ K , where K := λ (1 − ∆) G ( δ ) − λ ∆ . Now pick < k < n . If x k ≤ ∆ and x k − ≥ δ then f k = λ k − x k − (1 − x k ) − λ k x k (1 − x k +1 ) ≥ K k , where K k := λ k − δ (1 − ∆) − λ k ∆ . If x n ≤ ∆ and x n − ≥ δ for ≤ i ≤ n − then f n = λ n − x n − (1 − x n ) − λ n x n ≥ K n − . Finally, if x n +1 ≤ ∆ and x n ≥ δ then f n +1 = λ n x n − λ (1 − x ) G ( x n +1 ) ≥ K n +1 , where K n +1 := λ n δ − λ G (∆) . Thus, (26) holds for K := min { K , . . . , K n +1 } , and clearly K > forall ∆ > sufficiency small. Thus, (29) satisfies (CBR) .To show that (29) also satisfies (27), note that for all k ∈ { , . . . , n } and all x ∈ Ω , ˙ x k ≥ − λ k x k , andthat by the properties of G ( · ) , ˙ x n +1 ≥ − sx n +1 for all x n +1 > sufficiently small. Thus, (29) satisfies allthe conditions in Lemma 1 and this implies that for any τ > there exists ε = ε ( τ / > such that x i ( t, a ) ≥ ε, for all i ∈ { , . . . , n + 1 } and all t ≥ τ / . (30)This proves the first part of Proposition 1. To complete the proof, define y i ( t ) := 1 − x n +1 − i ( t ) , i = 1 , . . . , n, (31)and y n +1 ( t ) := x n +1 ( t ) . Then ˙ y = λ n (1 − y ) − λ n − y (1 − y ) , ˙ y = λ n − y (1 − y ) − λ n − y (1 − y ) , ... ˙ y n − = λ y n − (1 − y n − ) − λ y n − (1 − y n ) , ˙ y n = λ y n − (1 − y n ) − λ G ( y n +1 ) y n , ˙ y n +1 = λ n (1 − y ) − λ G ( y n +1 ) y n . The first n equations here are an RFM with a time-varying exit rate λ G ( y n +1 ( t )) . We already knowthat y n +1 ( t ) ≥ ε for all t ≥ τ / and the results in [39] imply that there exists ε = ε ( τ ) > such that y i ( t, a ) ≥ ε , for all i ∈ { , . . . , n } and all t ≥ τ. Combining this with (31) completes the proof of Prop. 1. (cid:3)
ProofofProposition2. The Jacobian matrix J of the RFMNP has the form J = J , ... J ,m +1 J , ... J ,m +1 ... ... ... ... ... J m,m J m,m +1 J m +1 , J m +1 , ... J m +1 ,m J m +1 ,m +1 ∈ R d × d . (32) Here J i,i ∈ R n i × n i , i = 1 , . . . , m , is the Jacobian of RFMIO i given by − λ i G i ( z ) − λ i (1 − x i ) λ i x i λ i (1 − x i ) − λ i x i − λ i (1 − x i ) λ i x i λ i (1 − x i ) − λ i x i − λ i (1 − x i ) 0 0 ... − λ in i − x in i − − λ in i − (1 − x in i ) λ in i − x in i − λ in i − (1 − x in i ) − λ in i − x in i − − λ in i , and the other blocks are J m +1 ,i = [ λ i G i ( z ) 0 ... λ in i ] ∈ R × n i ,J i,m +1 = [ λ i (1 − x i ) G ′ i ( z ) 0 ... ′ ∈ R n i × ,i = 1 , . . . , m , and J m +1 ,m +1 = − m X j =1 λ j (1 − x j ) G ′ j ( z ) ∈ R . Since x i ∈ [0 , and λ ij > , every off-diagonal entry of J i,i is non-negative. Since z ≥ , everyentry of J m +1 ,j , J i,m +1 , i = 1 , . . . , m , is also nonnegative. We conclude that every non-diagonal entryof J is non-negative, and this implies (11) (see [59]). To prove (12), note that for any point in Int(Ω) (i.e., x j ∈ (0 , , j = 1 , . . . , d − , and z > ) every entry on the super- and sub-diagonal of J i,i is strictlypositive. Also, G i ( z ) > , i = 1 , . . . , m . This implies that the matrix J in (32) is irreducible. Combiningthis with Prop. 1 completes the proof. (cid:3) Proof of Theorem 1. Since the RFMNP is a cooperative irreducible system on
Int(Ω) with a non-trivialfirst integral, Thm. 1 follows from combining Prop. 1 with the results in [46] (see also [51], [45] and [33]for related ideas). (cid:3)
ProofofProposition3. Recall that the matrix measure µ ( · ) : R n × n → R induced by the ℓ norm is µ ( A ) = max { c ( A ) , . . . , c n ( A ) } , where c i ( A ) := a ii + P k = i | a ki | , i.e. the sum of entries in column i of A , with the off-diagonal entriestaken with absolute value [68]. For the Jacobian of the RFMNP (32), we have c i ( J ( x )) = 0 for all i and all x ∈ Ω , so µ ( J ( x )) = 0 . Now (16) follows from standard results in contraction theory (see,e.g., [54]). (cid:3) Proof of Theorem 2. Write the PRFMNP as ˙ x = f ( t, x ) . Then f ( t, y ) = f ( t + T, y ) for all t and y .Furthermore, H in (9) is a non trivial first integral of the dynamics. Now Theorem 2 follows from theresults in [61] (see also [30]). The fact that φ s ∈ Int(Ω) follows from Proposition 1. (cid:3)
Proof of Theorem 3. To simplify the presentation, we prove the theorem for the case m = 2 and changesome of the notation. The proof in the general case is very similar. Let n [ ℓ ] denote the dimension ofthe first [second] RFMIO, and let λ i [ ζ i ] denote the rates in the first [second] RFMIO. We denote thestate-variables of RFMIO by x i , i = 1 , . . . , n , and those of the second RFMIO by y i , i = 1 , . . . , ℓ .Let e i , i = 1 , . . . , n , [ v j , j = 1 , . . . , ℓ ] denote the equilibrium point of the first [second] RFMIO. Thenwe need to prove that ˜ e i < , and ˜ e j > , for all j ∈ { i + 1 , ..., n } , (33) sign(˜ v ) = · · · = sign(˜ v ℓ ) = sign(˜ e z ) , (34) where ˜ v j := ¯ v j − v j . At steady state, the RFMNP equations yield: λ G ( e z )(1 − e ) = λ e (1 − e ) , = λ e (1 − e ) , = λ e (1 − e ) , ... = λ n − e n − (1 − e n )= λ n e n , (35) ζ G ( e z )(1 − v ) = ζ v (1 − v ) , = ζ v (1 − v ) , = ζ v (1 − v ) , ... = ζ ℓ − v ℓ − (1 − v ℓ )= ζ ℓ v ℓ , (36) λ n e n + ζ ℓ v ℓ = λ G ( e z )(1 − e ) + ζ G ( e z )(1 − v ) . (37)Also, since H is a first integral n X i =1 e i + ℓ X j =1 v j + e z = H (0) . (38)Pick i ∈ { , . . . , n − } . Consider a new RFMNP obtained by modifying λ i in RFMIO to ¯ λ i ,with ¯ λ i > λ i . Then the equations for the modified equilibrium point are: λ G (¯ e z )(1 − ¯ e ) = λ ¯ e (1 − ¯ e ) , = λ ¯ e (1 − ¯ e ) , ... = ¯ λ i ¯ e i (1 − ¯ e i +1 ) , ... = λ n − ¯ e n − (1 − ¯ e n )= λ n ¯ e n , (39) ζ G (¯ e z )(1 − ¯ v ) = ζ ¯ v (1 − ¯ v ) , = ζ ¯ v (1 − ¯ v ) , ... = ζ ℓ − ¯ v ℓ − (1 − ¯ v ℓ )= ζ ℓ ¯ v ℓ , (40) λ n ¯ e n + ζ ℓ ¯ v ℓ = λ G (¯ e z )(1 − ¯ e ) + ζ G (¯ e z )(1 − ¯ v ) . (41)Since the initial condition remains the same, n X i =1 e i + ℓ X j =1 v j + e z = n X i =1 ¯ e i + ℓ X j =1 ¯ v j + ¯ e z = H (0) , so n X i =1 ˜ e i + ℓ X j =1 ˜ v j + ˜ e z = 0 . (42)The last equality in (36) yields ζ ℓ − v ℓ − = ζ ℓ v ℓ − v ℓ . The right-hand side here is increasing in v ℓ , and the left-hand side is increasing in v ℓ − (recall that ζ ℓ and ζ ℓ − are the same for both the original and the modified RFMNP), so a change in λ i must leadto sign(˜ v ℓ − ) = sign(˜ v ℓ ) . Using (36) again yields ζ ℓ − v ℓ − = ζ ℓ v ℓ − v ℓ − , so sign(˜ v ℓ − ) = sign(˜ v ℓ − ) = sign(˜ v ℓ ) . Continuing in this way yields sign(˜ v ) = sign(˜ v ) = · · · = sign(˜ v ℓ ) . (43)By (36), G ( e z ) = ζ ℓ v ℓ ζ (1 − v ) . Since G i ( p ) is strictly increasing in p , combining this with (43) implies that sign(˜ e z ) = sign(˜ v ) =sign(˜ v ) = · · · = sign(˜ v ℓ ) . This proves (34). To prove (33), note that arguing as above using (39) yields sign(˜ e n ) = sign(˜ e n − ) = · · · = sign(˜ e i +1 ) . Seeking a contradiction, assume that sign(˜ e i ) = sign(˜ e n ) . By (39), e i − = λ n e n λ i − (1 − e i ) so sign(˜ e i − ) = sign(˜ e n ) , and continuing in this fashion yields sign(˜ e n ) = sign(˜ e n − ) = · · · = sign(˜ e ) . (44)By (35), G ( e z ) = λ n e n λ (1 − e ) , and combining this with (44) implies that sign(˜ e ) = ... = sign(˜ e n ) = sign(˜ e z ) . Since also sign(˜ v ) = ... =sign(˜ v ℓ ) = sign(˜ e z ) , it follows that all the differences have the same sign, and this contradicts (42). Weconclude that if ˜ e i = 0 then sign(˜ e i ) = sign(˜ e i +1 ) = sign(˜ e i +2 ) = · · · = sign(˜ e n ) . To complete the proofof (33), we need to show that ˜ e i < . Seeking a contradiction, assume that ˜ e i ≥ , so ˜ e i +1 ≤ , . . . , ˜ e n ≤ .Thus, ¯ e i ≥ e i , ¯ e i +1 ≤ e i +1 , . . . , ¯ e n ≤ e n . By (35), λ i = λ n e n e i (1 − e i +1 ) , so ¯ λ i − λ i λ n = ¯ e n ¯ e i (1 − ¯ e i +1 ) − e n e i (1 − e i +1 ) ≤ . This contradiction completes the proof for the case where i ∈ { , . . . , n − } . The proof for the case i = 0 and i = n is similar. (cid:3) Proof of Proposition 4. In [45], Mierczynski considered an irreducible cooperative dynamical sys-tem, ˙ x = f ( x ) , that admits a non trivial first integral H ( x ) with a positive gradient. Let S x := { v ∈ R n : ∇ H T ( x ) v = 0 } , and consider the extended system ˙ x = f ( x ) , ˙ δx = J ( x ) δx, where J := ∂∂x f is the Jacobian, with initial condition x (0) = x , δx (0) = δx ∈ S x \ { } . Mierczynskishows that there exists a norm, that depends on x , such that | δx ( t ) | x ( t ) < | δx | x , for all t > . (For a general treatment on using Lyapunov-Finsler functions in contraction theory, see [18].) At theunique equilibrium point e this yields | exp( J ( e ) t ) δx | e < | δx | e , for all t > . This implies that the matrix obtained by restricting J ( e ) to the integral manifold is Hurwitz and thusnonsingular. Invoking the implicit function theorem implies that the mapping from λ to e can be identifiedwith a C function. (cid:3) R EFERENCES [1] D. Adams, B. Schmittmann, and R. Zia, “Far-from-equilibrium transport with constrained resources,”
J. Stat. Mech. , p. P06009, 2008.[2] U. Ala, F. A. Karreth, C. Bosia, A. Pagnani, R. Taulli, V. Leopold, Y. Tay, P. Provero, R. Zecchina, and P. P. Pandolfi, “Integratedtranscriptional and competitive endogenous RNA networks are cross-regulated in permissive molecular environments,”
Proceedings ofthe National Academy of Sciences , vol. 110, no. 18, pp. 7154–9, 2013.[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] D. Angeli and E. D. Sontag, “Monotone control systems,”
IEEE Trans. Automat. Control , vol. 48, pp. 1684–1698, 2003.[5] T. Ben-Yehezkel, S. Atar, H. Zur, A. Diament, E. Goz, T. Marx, R. Cohen, A. Dana, A. Feldman, E. Shapiro, and T. Tuller, “Rationallydesigned, heterologous S. cerevisiae transcripts expose novel expression determinants,”
RNA Biol. , 2015.[6] R. A. Blythe and M. R. Evans, “Nonequilibrium steady states of matrix-product form: a solver’s guide,”
J. Phys. A: Math. Theor. ,vol. 40, no. 46, pp. R333–R441, 2007.[7] C. A. Brackley, M. C. Romano, and M. Thiel, “The dynamics of supply and demand in mRNA translation,”
PLoS Comput. Biol. ,vol. 7, no. 10, p. e1002203, 10 2011.[8] H. Bremer and P. P. Dennis,
Modulation of chemical composition and other parameters of the cell by growth rate , 2nd ed. Escherichiacoli and Salmonella typhimurium: Cellular and Molecular Biology, 1996, ch. 97, p. 1559.[9] R. C. Brewster, F. M. Weinert, H. G. Garcia, D. Song, M. Rydenfelt, and R. Phillips, “The transcription factor titration effect dictateslevel of gene expression,”
Cell , vol. 156, no. 6, pp. 1312–1323, 2014.[10] L. Carmel and E. Koonin, “A universal nonmonotonic relationship between gene compactness and expression levels in multicellulareukaryotes,”
Genome Biol Evol. , vol. 1, pp. 382–90, 2009.[11] C. Castillo-Davis, S. Mekhedov, D. Hartl, E. Koonin, and F. Kondrashov, “Selection for short introns in highly expressed genes.”
NatGenet. , vol. 31, no. 4, pp. 415–8, 2002.[12] L. J. Cook and R. K. P. Zia, “Feedback and fluctuations in a totally asymmetric simple exclusion process with finite resources,”
J. Stat.Mech. , p. P02012, 2009.[13] L. J. Cook and R. K. P. Zia, “Competition for finite resources,”
J. Stat. Mech. , p. P05008, 2012.[14] A. Dana and T. Tuller, “Efficient manipulations of synonymous mutations for controlling translation rate–an analytical approach.”
J.Comput. Biol. , vol. 19, pp. 200–231, 2012.[15] A. Dana and T. Tuller, “The effect of trna levels on decoding times of mrna codons,”
Nucleic Acids Res. , vol. 42, no. 14, pp. 9171–81,2014.[16] D. Del Vecchio and R. M. Murray,
Biomolecular Feedback Systems . Princeton University Press, 2014.[17] E. Eisenberg and E. Levanon, “Human housekeeping genes are compact,”
Trends Genet. , vol. 19, no. 7, pp. 362–5, 2003.[18] F. Forni and R. Sepulchre, “A differential Lyapunov framework for contraction analysis,”
IEEE Trans. Automat. Control , vol. 59, no. 3,pp. 614–628, 2014.[19] M. Frenkel-Morgenstern, T. Danon, T. Christian, T. Igarashi, L. Cohen, Y. M. Hou, and L. J. Jensen, “Genes adopt non-optimal codonusage to generate cell cycle-dependent oscillations in protein levels,”
Mol. Syst. Biol. , vol. 8, p. 572, 2012.[20] T. Galla, “Optimizing evacuation flow in a two-channel exclusion process,”
J. Stat. Mech. , p. P09004, 2011.[21] P. Greulich, L. Ciandrini, R. J. Allen, and M. C. Romano, “Mixed population of competing totally asymmetric simple exclusionprocesses with a shared reservoir of particles,”
Phys. Rev. E , vol. 85, p. 011142, 2012.[22] A. Gyorgy and D. Del Vecchio, “Limitations and trade-offs in gene expression due to competition for shared cellular resources,” in
Proc. 53rd IEEE Conf. on Decision and Control , Dec. 2014, pp. 5431–5436. [23] A. Gyorgy, J. I. Jimenez, J. Yazbek, H. Huang, H. Chung, R. Weiss, and D. Del Vecchio, “Isocost lines describe the cellular economyof genetic circuits,” Biophysical J. , 2015, To appear.[24] M. Ha and M. den Nijs, “Macroscopic car condensation in a parking garage,”
Phys. Rev. E , vol. 66, p. 036118, 2002.[25] A. Hamadeh and D. Del Vecchio, “Mitigation of resource competition in synthetic genetic circuits through feedback regulation,” in
Proc. 53rd IEEE Conf. on Decision and Control , Dec. 2014, pp. 3829–3834.[26] Q. Hui and W. M. Haddad, “Distributed nonlinear control algorithms for network consensus,”
Automatica , vol. 44, no. 9, pp. 2375–2381,2008.[27] N. T. Ingolia, S. Ghaemmaghami, J. R. Newman, and J. S. Weissman, “Genome-wide analysis in vivo of translation with nucleotideresolution using ribosome profiling,”
Science , vol. 324, no. 5924, pp. 218–23, 2009.[28] N. Jacques and M. Dreyfus, “Translation initiation in Escherichia coli: old and new questions,”
Mol Microbiol. , vol. 4, no. 7, pp.1063–7, 1990.[29] M. Jens and N. Rajewsky, “Competition between target sites of regulators shapes post-transcriptional gene regulation,”
Nat Rev Genet. ,vol. 16, no. 2, pp. 113–26, 2015.[30] J. Ji-Fa, “Periodic monotone systems with an invariant function,”
SIAM J. Math. Anal. , vol. 27, pp. 1738–1744, 1996.[31] G. Kopeina, Z. A. Afonina, K. Gromova, V. Shirokov, V. Vasiliev, and A. Spirin, “Step-wise formation of eukaryotic double-rowpolyribosomes and circular translation of polysomal mRNA,”
Nucleic Acids Res. , vol. 36, no. 8, pp. 2476–2488, 2008.[32] G. Kudla, A. W. Murray, D. Tollervey, and J. B. Plotkin, “Coding-sequence determinants of gene expression in Escherichia coli,”
Science , vol. 324, pp. 255–258, 2009.[33] P. D. Leenheer, D. Angeli, and E. D. Sontag, “Monotone chemical reaction networks,”
J. Mathematical Chemistry , vol. 41, pp. 295–314,2007.[34] G.-W. Li, E. Oh, and J. Weissman, “The anti-Shine-Dalgarno sequence drives translational pausing and codon choice in bacteria,”
Nature , vol. 484, no. 7395, pp. 538–41, 2012.[35] S. Li, L. Feng, and D. Niu, “Selection for the miniaturization of highly expressed genes,”
Biochem Biophys Res Commun. , vol. 360,no. 3, pp. 586–92, 2007.[36] W. Lohmiller and J.-J. E. Slotine, “On contraction analysis for non-linear systems,”
Automatica , vol. 34, pp. 683–696, 1998.[37] J. Lu and C. Deutsch, “Electrostatics in the ribosomal tunnel modulate chain elongation rates,”
J. Mol. Biol. , vol. 384, p. 736, 2008.[38] C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin, “Kinetics of biopolymerization on nucleic acid templates,”
Biopolymers , vol. 6, pp.1–25, 1968.[39] 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.[40] M. Margaliot, E. Sontag, and T. Tuller, “Contraction after small transients,” Submitted. [Online]. Available: arXiv:1506.06613[41] 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.[42] M. Margaliot and T. Tuller, “Stability analysis of the ribosome flow model,”
IEEE/ACM Trans. Computational Biology andBioinformatics , vol. 9, pp. 1545–1552, 2012.[43] Margaliot, M. and Tuller, T., “Ribosome flow model with positive feedback,”
J. Royal Society Interface , vol. 10, p. 20130267, 2013.[44] W. Mather, J. Hasty, L. Tsimring, and R. Williams, “Translational cross talk in gene networks,”
Biophys J. , vol. 104, no. 11, pp.2564–72, 2013.[45] J. Mierczynski, “A class of strongly cooperative systems without compactness,”
Colloq. Math. , vol. 62, pp. 43–47, 1991.[46] J. Mierczynski, “Cooperative irreducible systems of ordinary differential equations with first integral,”
ArXiv e-prints , 2012. [Online].Available: http://arxiv.org/abs/1208.4697[47] S. Pechmann and J. Frydman, “Evolutionary conservation of codon optimality reveals hidden signatures of cotranslational folding,”
Nat Struct Mol Biol. , vol. 20, no. 2, pp. 237–43, 2013.[48] 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.[49] G. Poker, M. Margaliot, and T. Tuller, “Sensitivity of mRNA translation,”
Sci. Rep. , vol. 5, p. 12795, 2015.[50] Y. Qian and D. Del Vecchio, “Effective interaction graphs arising from resource limitations in gene networks,” in
Proc. AmericanControl Conference (ACC2015) , Chicago, IL, 2015.[51] A. Raveh, Y. Zarai, M. Margaliot, and T. Tuller, “Ribosome flow model on a ring,”
IEEE/ACM Trans. Computational Biology andBioinformatics , 2015, To appear. [Online]. Available: http://ieeexplore.ieee.org/xpl/articleDetails.jsp?arnumber=7078877[52] 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.[53] J. D. Richter and L. D. Smith, “Differential capacity for translation and lack of competition between mRNAs that segregate to freeand membrane-bound polysomes,”
Cell , vol. 27, pp. 183–191, 1981.[54] G. Russo, M. di Bernardo, and E. D. Sontag, “Global entrainment of transcriptional systems to periodic inputs,”
PLOS ComputationalBiology , vol. 6, p. e1000739, 2010.[55] R. Sabi and T. Tuller, “A comparative genomics study on the effect of individual amino acids on ribosome stalling,”
BMC genomics ,2015.[56] B. Schwanhausser, D. Busse, N. Li, G. Dittmar, J. Schuchhardt, J. Wolf, W. Chen, and M. Selbach, “Global quantification of mammaliangene expression control,”
Nature , vol. 473, no. 7347, pp. 337–42, 2011.[57] P. M. Sharp, L. R. Emery, and K. Zeng, “Forces that influence the evolution of codon bias,”
Philos. Trans. R. Soc. Lond. B , vol. 365,no. 1544, pp. 1203–12, 2010.[58] 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. [59] 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.[60] P. Spellman, G. Sherlock, M. Zhang, V. Iyer, K. Anders, M. Eisen, P. Brown, D. Botstein, and B. Futcher, “Comprehensive identificationof cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization,”
Mol Biol Cell. , vol. 9, no. 12, pp.3273–97, 1998.[61] B. Tang, Y. Kuang, and H. Smith, “Strictly nonautonomous cooperative system with a first integral,”
SIAM J. Math. Anal. , vol. 24, pp.1331–1339, 1993.[62] G. Tripathy and M. Barma, “Driven lattice gases with quenched disorder: exact results and different macroscopic regimes,”
Phys. Rev.E , vol. 58, pp. 1911–1926, 1998.[63] T. Tuller, A. Carmi, K. Vestsigian, S. Navon, Y. Dorfan, J. Zaborske, T. Pan, O. Dahan, I. Furman, and Y. Pilpel, “An evolutionarilyconserved mechanism for controlling the efficiency of protein translation,”
Cell , vol. 141, no. 2, pp. 344–54, 2010.[64] T. Tuller, E. Ruppin, and M. Kupiec, “Properties of untranslated regions of the S. cerevisiae genome,”
BMC Genomics , vol. 10, p. 329,2009.[65] T. Tuller, I. Veksler, N. Gazit, M. Kupiec, E. Ruppin, and M. Ziv, “Composite effects of the coding sequences determinants on thespeed and density of ribosomes,”
Genome Biol. , vol. 12, no. 11, p. R110, 2011.[66] T. Tuller, Y. Y. Waldman, M. Kupiec, and E. Ruppin, “Translation efficiency is determined by both codon bias and folding energy,”
Proceedings of the National Academy of Sciences , vol. 107, no. 8, pp. 3645–50, 2010.[67] 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.[68] M. Vidyasagar,
Nonlinear Systems Analysis . Englewood Cliffs, NJ: Prentice Hall, 1978.[69] J. Vind, M. A. Sorensen, M. D. Rasmussen, and S. Pedersen, “Synthesis of proteins in Escherichia coli is limited by the concentrationof free ribosomes: expression from reporter genes does not always reflect functional mRNA levels,”
J. Mol. Biol. , vol. 231, pp. 678–88,1993.[70] T. von der Haar, “Mathematical and computational modelling of ribosomal movement and protein synthesis: an overview,”
Computationaland Structural Biotechnology J. , vol. 1, no. 1, pp. 1–7, 2012.[71] J. Warner, “The economics of ribosome biosynthesis in yeast,”
Trends Biochem Sci. , vol. 24, no. 11, pp. 437–40, 1999.[72] M. Welch, S. Govindarajan, J. E. Ness, A. Villalobos, A. Gurney, J. Minshull, and C. Gustafsson, “Design parameters to controlsynthetic gene expression in Escherichia coli,”
PLoS One , vol. 4, no. 9, p. e7002, 2009.[73] Y. Zarai, M. Margaliot, and T. Tuller, “Maximizing protein translation rate in the ribosome flow model: the homogeneous case,”
IEEE/ACM Trans. Computational Biology and Bioinformatics , vol. 11, pp. 1184–1195, 2014.[74] D. Zenklusen, D. R. Larson, and R. H. Singer, “Single-RNA counting reveals alternative modes of gene expression in yeast,”