key: cord-0304561-pw91t3v1 authors: Jain, Aditi; Margaliot, Michael; Gupta, Arvind Kumar title: Large-scale mRNA translation and the intricate effects of competition for the finite pool of ribosomes date: 2021-09-16 journal: bioRxiv DOI: 10.1101/2021.09.15.460428 sha: 419a0eaacb079d871faa8bb668e42598f95d8447 doc_id: 304561 cord_uid: pw91t3v1 We present a new theoretical framework for large-scale mRNA translation using a network of models called the ribosome flow model with Langmuir kinetics (RFMLK), interconnected via a pool of free ribosomes. The input to each RFMLK depends on the pool density, and it affects the initiation rate and the internal ribosome entry rates at each site along each RFMLK. Ribosomes that detach from an RFMLK due to termination or premature drop-off are fed back into the pool. We prove that the network always converges to a steady-state, and study its sensitivity to variations in the parameters. For example, we show that if the drop-off rate at some site in some RFMLK is increased then the pool density increases and consequently the steady-state production rate in all the other RFMLKs increases. Surprisingly, we also show that modifying a parameter of a certain RFMLK can lead to arbitrary effects on the densities along the modified RFMLK, depending on the parameters in the entire network. We conclude that the competition for shared resources generates an indirect and intricate web of mutual effects between the mRNA molecules, that must be accounted for in any analysis of translation. Gene expression is a crucial cellular process that produces desired proteins by decoding the genetic code in the DNA [1] . In Eukaryotes, gene expression includes several processes. During transcription, the protein-coding information encoded in a DNA is copied into an mRNA molecule. During translation, the information inscribed in the mRNA is translated to chains of amino acids destined to become proteins [2, 3] . The information encoded in the DNA and mRNA is decoded by complex molecular machines called RNA polymerases and ribosomes, respectively. The dynamical flow of these machines along the DNA or mRNA "tracks" plays an important role in gene expression and, in particular, in regulating protein production. Translation is a fundamental process that takes place in all living cells, from bacteria to humans. Many mRNA molecules are translated in parallel with many ribosomes decoding every mRNA molecule [4, 5] . This implies that the mRNA molecules effectively "compete" for the finite resources in the cell, like tRNA molecules and free ribosomes [6, 7] . The competition for ribosomes may explain important dynamical properties of translation, that are difficult to understand when considering the translation of a single, isolated mRNA molecule. For example, it is known that stalling ribosomes may detach from the mRNA before completing the translation process [8, 9] . This is somewhat surprising, as a ribosome that drops-off from the mRNA before reaching the stop codon fails to complete the synthesis of a full-length protein, and releases a truncated protein, whose accumulation could be detrimental to the cell [10] . However, in the context of competition for free ribosomes, premature drop-off may have a positive effect: it allows stalled ribosomes to join the pool of free ribosomes that can initiate translation in other mRNA molecules. Thus, modeling translation as a network of interconnected processes and taking into account competition for shared resources is important for gaining a deeper understanding of fundamental principles in cellular biophysics. Typically, translation is initiated by ribosome scanning from the 5 end of the capped mRNA. However, some mRNAs include internal ribosome entry sites (IRESs), that allow for translation initiation in a cap-independent manner. IRESs, first discovered in poliovirus, are common in RNA viruses and allow viral translation even when host translation is inhibited for some reason [11, 14] . Cellular growth regulatory genes, and genes transcribed in response to stress also contain IRES elements [12] . Also, synthetic biologists often insert IRES sequences into their vectors to allow expressing two or more genes from a single vector [15] . We believe that the effect of IRESs on translation should also take into account the competition for shared resources, like ribosomes. For example, a recent study [13] shows that the non-structural protein 1 (Nsp1), produced by the SARS-CoV-2 virus, binds to the human 40S subunit in ribosomal complexes, and thus interferes with mRNA binding. Since translation of the viral mRNA is more efficient than that of cellular mRNAs with 5 UTRs, the net effect is hijacking the cellular translation machinery by the virus. Many mathematical and computational models for translation on a single mRNA molecule have been developed [16] . One such model is the Asymmetric Simple Exclusion Process (ASEP), and its variants motivated by biological findings [17] [18] [19] [20] [21] [22] [23] . This model includes a 1D chain of sites, and particles that hop according to a stochastic mechanism from a site to a neighboring conditions corresponding to an equal total number of ribosomes in the network converge to the same equilibrium point. In other words, the network "forgets" the exact initial condition, except for the total initial density of ribosomes. This qualitative behaviour holds for any feasible set of parameters covering many possible biophysical conditions. We also show that if all the transition rates vary in a periodic manner, with a common period T , then every solution of the RFMLKN converges to a periodic solution with period T . This implies in particular that the protein production rate entrains to periodic excitations in the translational machinery. We use the RFMLKN to analyze quantitatively and qualitatively important questions such as the effect of ribosome drop-off/attachment from/to a site in one mRNA on the steady-state production rate of all the other mRNAs in the network. The analysis highlights how the competition for shared resources generates an indirect and intricate web of mutual effects between the mRNA molecules in the cell. These effects cannot be analyzed using models of translation on a single, isolated mRNA molecule. At each elongation step, every sequence of three consecutive nucleotides in the mRNA, called a codon, is decoded into an amino-acid, and this process continues until the ribosome reaches the 3 end of the mRNA [2] . The codon decoding rates may vary among different mRNAs and depend on many transcript features [40] . Several ribosomes may scan the same mRNA molecule in parallel, but a ribosome cannot overtake another ribosome in front of it, thus obeying the simple exclusion principle. Ribosomes may detach from the mRNA molecule before reaching the stop codon due to several reasons like ribosome "traffic jams", the presence of a premature stop codon, ribosome-ribosome interactions due to depletion of aminoacyl tRNA or amino-acid misincorporation, etc. [41] . The limited availability of free ribosomes induces indirect coupling due to competition between mRNA molecules. We prove that for a given set of elongation, drop-off and attachment rates, and a total number of ribosomes in the network, the RFMLKN admits a unique steady-state i.e. the ribosomal density profiles on all the mRNAs and in the pool converge to a fixed value, as time goes to infinity. This raises the important question of how does the steady-state changes if we modify any parameter in the model, e.g. the rate of ribosome drop-off in one site of a specific mRNA. Our analysis shows that following an increase (decrease) in a drop-off rate in any mRNA molecule, the steady-state ribosome density profile in all the other mRNA molecules increases (decreases). The intuitive explanation for this is as follows: increasing the drop-off rate leads to releasing more ribosomes to the pool of free ribosomes and this increases the initiation rate as well as the attachment rate in all the other mRNA molecules leading to an increase in their ribosome density profile. We also prove the "dual" result, namely, that increasing (decreasing) an attachment rate in a specific mRNA decreases (increases) the steady-state ribosome densities in all the other mRNA molecules. However, and perhaps surprisingly, we show that it is very difficult to predict the effect of a variation in one of the rates on the mRNA that is modified, as the effect will depend on the entire network. For example, increasing the attachment rate in one site of a specific mRNA may deplete the pool and thus decrease the effective attachment rates in other sites along the modified mRNA, leading to an unexpected decrease in the density along this mRNA. These results highlight the indirect effects of competition for resources, and also the importance of taking competition into account when "tinkering" with the bio-physical features of a single mRNA molecule e.g. by replacing codons by synonymous codons. Our simulations suggest another interesting implication of ribosome drop-off and/or ribosome attachment (e.g. in IRESs). It seems that these phenomena increase the amount of indirect "communication" between the mRNA molecules, through the pool, and thus lead to a higher level of synchronization between the production rates in the mRNAs. This suggests another possible advantage of ribosome drop-off and attachment as tools for regulating the protein production from different copies of the same mRNA. Our model is a network of interconnected ribosome flow models with Langmuir kinetics (RFMLK). We begin by reviewing the RFMLK and then describe the network of interconnected RFMLKs. The RFMLK is a deterministic, non-linear, continuous-time compartmental model for modeling ribosome flow along a single mRNA molecule and is a coarse-grained mean field approximation of TASEP with Langmuir kinetics and open boundary conditions [28, 35] . Figure 1 : The RFMLK models unidirectional flow along a chain of n sites. The density at site i at time t is represented by x i (t). The parameter λ i > 0 controls the transition rate from site i to site i + 1, with λ 0 and λ n controlling the initiation and termination rates, respectively. The parameter . . . Defining x 0 (t) ≡ 1 and x n+1 (t) ≡ 0 allows to write these equations in the more succinct form: The term λ i−1 x i−1 (1 − x i ) represents the flow of the particles from site i − 1 to site i. This increases with the density at site i − 1 and decreases as site i becomes fuller. This corresponds to a "soft" version of the simple exclusion principle in ASEP. Similarly, the term represents the flow from site i to i + 1. The term β i (1 − x i ) represents the attachment of particles from the environment to site i, whereas α i x i represents the drop-off of particles from site i to the environment. By setting some of the α i s and β i s to positive values and the others to zero, it is possible to model drop-off and attachment at specific sites. If α i = β i = 0 for all i then the model reduces to the RFM. The topology of the RFMLK is depicted in Figure 1 . To build a network, we use an RFMLK with an input and output given bẏ . . The time-varying control u(t) multiplies the term representing the entry rate into site 1, and also the attachment rates in all the sites. We assume that u(t) ≥ 0 for any time t. A larger value of u(t) corresponds for example to a larger density of free ribosomes in the vicinity of the mRNA at time t, and consequently it increases the effective initiation rate in the first site and the attachment rates in all the sites. The output y(t) is the total exit rate of ribosomes from the RFMLK to the environment at time t. Note that (3) is a nonlinear model, as it includes both products of state-variables and products of state-variables and the control input. Example 1 Consider the case n = 1. In this special case, (3) becomes the linear systeṁ (0, 1) for all t > 0, and that the limit e 1 := lim t→∞ x 1 (t) exists, and satisfies In particular, if v = 0 then e 1 = 0, and if v → ∞ then e 1 → 1. The first case corresponds to no ribosomes in the vicinity of the mRNA, so the single site empties. The second case corresponds to an infinite density of ribosomes, so the site fills up completely. Note also that e 1 is an increasing function of λ 0 , β 1 , v, and a decreasing function of λ 1 , α 1 . The next subsection introduces the RFMLKN. To model competition for the finite pool of ribosomes in the cell, we consider a set of m RFMLKs with input and output, representing m different mRNA molecules in the cell, connected via a pool of free ribosomes. The ith RFMLK has length n i , input function u i , output y i , and rates λ i The dynamics of the m RFMLKs is written aṡ These RFMLKs are interconnected through a pool of free ribosomes, i.e. ribosomes that are not attached to any mRNA. We use the scalar function z(t) ≥ 0 to denote the density of ribosomes in the pool at time t. The pool feeds the initiation locations as well as the sites in the mRNAs where attachment takes place. Mathematically, this implies that u i (t) = G i (z(t)), i = 1, 2, . . . , m. We assume that every function G i (·) : R + → R + satisfies the following properties: 2. G i (·) is continuously differentiable with G i (z) > 0 for all z ≥ 0; and 3. There exists c > 0 such that G i (z) ≤ cz for all z > 0 sufficiently small. The first property implies that if the pool is empty then no ribosomes can exit the pool; the second implies that as the number of ribosomes in the pool increases, more ribosomes exit the pool; and the third is a technical condition that is needed for the proof given later on of persistence in the RFMLKN. Functions that satisfy these properties include, for example, the linear function G(z) = az, with a > 0, and the bounded function G(z) = a tanh(bz), with a, b > 0 [20, 42] . The dynamics of the ith RFMLK in the network is thus given by: Figure 3 : Each mRNA is described by an RFMLK with input and output. The output of each RFMLK is fed into the pool, and the pool feeds the initiation and attachment rates in all the RFMLKs. The output of each RFMLK is fed into the pool. Hence, the pool dynamics is given by: In other words, all the ribosomes that exit or drop-off the mRNAs feed the pool, and the pool feeds the initiation and attachment sites in all the mRNAs. Summarizing, the RFMLKN is a dynamical system with d := 1+ m i=1 n i state variables, and dynamics described by equations (4), (5) and (6) (see Figure 3 ). Eq. (6) and our assumptions on the functions G i imply that if z(0) ≥ 0 then z(t) ≥ 0 for all t ≥ 0, i.e. the pool density is always non-negative. Let This is the total number of ribosomes in the system at time t. An important property of the RFMLKN is that it is a closed system, so Thus, H is a first integral of the dynamics. The RFMLKN models the indirect coupling between the mRNA molecules induced by competition for the finite number of ribosomes in the system. For example, if there is a "traffic jam" of ribosomes on one of the mRNAs then the pool is depleted and thus the initiation and attachment rates to all the mRNAs decrease. We prove in Section 4 that all the state variables in the RFMLKN converge to a steady-state. Our first example demonstrates how the total number of ribosomes in the system affects the ribosomal densities in the mRNAs. Example 2 Consider an RFMLKN that includes a single RFMLK with dimension n 1 = 3, rates , and a pool with an output function G(z) = tanh(z). We simulated this system with the initial condition x 1 j = 0 for all j, and z(0) = c, so that the total number of ribosomes in the system is c, for various values of c. When c = 0 there are no ribosomes in the network and the steady-state values are all zero. As c increases, the number of ribosomes along the RFMLK increases. Since tanh(z) → 1 as z → ∞, the initiation and attachment rates converge to 1. Thus, ribosomal densities in RFMLK saturate to the values corresponding to initiation rate λ 1 0 = 1, and attachment rates β 1 i = 1. The remaining ribosomes accumulate in the pool (see Figure 4a ). Using a different output function, namely, G(z) = z, we see in Figure 4b that as c increases, all ribosomal densities tend to the maximal possible value 1, that is, the RFMLK completely fills up and the remaining ribosomes accumulate in the pool. Example 3 Consider an RFMLKN with m = 3 RFMLKs with dimensions n i = 3, i = 1, 2, 3, and parameters λ 1 We vary the parameter α, i.e. the ribosome drop-off rate from all the sites in the first RFMLK, in the range [0. 1, 3] . Figure 5 depicts the ASSD in each RFMLK as a function of α. It can be seen that as α increases, the ASSD in the first RFMLK decreases, whereas the ASSD in all the other RFMLKs increases. Indeed, as the drop-off rate from the first RFMLK increases, the density in the pool increases, and more ribosomes become available for translating the other mRNA molecules, thus increasing the ASSD in the other RFMLKs. From a biological perspective, this example corresponds to a situation when due to genetic errors or insufficient availability of charged tRNAs or frameshifting [43] , ribosomes start detaching before reaching the stop codon in an mRNA, resulting in truncated protein products. Our results explain why this may still be beneficial to the cell. The ribosome drop-off from one mRNA molecule increases the number of free ribosomes that are now available to translate other mRNAs which in turn increases the corresponding protein production rates. The next example considers the "dual" case of increasing the attachment rate in one of the mRNA molecules in the network. Consider an RFMLKN with m = 2 RFMLKs with dimensions n i = 10, i = 1, 2. The parameter values are λ 1 The initial condition is x i j = 0, for all i, j, and z(0) = 3.5. Figure 6 depicts the ASSD in each RFMLK as a function of β ranging from 0.1 to 3. It can be seen that as β increases, the ASSD in the first RFMLK increases and the ASSD in the other RFMLK decreases. This is due to the attachment of ribosomes at the first RFMLK leading to a depletion of ribosomes in the pool, and thus to a decrease in the ASSD in the second RFMLK. Note the relatively sharp decrease in the steady-state pool density as β increases. This is due to the fact that the number of sites is n i = 10, so a "traffic jam" in an RFMLK involves many stalled ribosomes along the RFMLK. (2) concomitantly slows down the cellular innate immune response [47] . IRESs have also been used as a biotechnological tool allowing the synthesis of several proteins of interest from one multicistronic mRNA [44] [45] [46] . In this context, the example above shows that the design of such tools should also take into consideration their effect on the pool of free ribosomes. The next example demonstrates the effect of modifying the length of one mRNA molecule in the network. Example 5 Consider an RFMLKN with m = 2 RFMLKs with dimensions n 1 = 5, and n 2 , respectively. The parameter values are λ i 0 = 1, λ i j = 1, α i j = 0.1, β i j = 0.01, for all i, j, and G i (z) = tanh(z), i = 1, 2. The initial condition is x i j = 0, for all i, j, and z(0) = 5. We simulated this network for various values of n 2 . As n 2 increases, there is a decrease in the ASSD in both RFMLKs and in the pool density (see Figure 7) . Indeed, increasing n 2 implies that ribosomes that bind to the second RFMLK remain on it for a longer period of time. This decreases the steady-state pool density and, consequently, the steady-state densities in all the RFMLKs. The next example again studies the effect of increasing n 2 and also compares the RFMLKN and the RFMNP. equal values. However, ρ 1 n 1 and ρ 2 n 2 are different. The reason for this may be that the non-zero attachment and detachment rates increase the indirect communication between the RFMLKs (through the pool) leading to better "synchronization". The next section rigorously analyzes the RFMLK and the RFMLKN using tools from systems and control theory and in particular the theory of cooperative dynamical systems [48] . We begin by analyzing the properties of the RFMLK with input and output described in (3), as these are the basic ingredients of the RFMLKN. Recall that x i (t) ∈ [0, 1] for all t, so the statespace of the RFMLK is C n := [0, 1] n . Let int(C n ) and ∂C n denote the interior and boundary of C n respectively. Let x(t, a) denote the solution of equation (3) at time t ≥ 0 for the initial condition a ∈ C n . For the sake of readability, all the proofs are placed in the Appendix. If u(t) ≡ 0 then no ribosomes enter the RFMLK and then it is clear that x(t) will converge to zero, that is, the density of ribosomes at each site will go to zero. The next result shows that for any input that is bounded below by a positive number all the state-variables remain bounded away from zero and also bounded away from one. For any τ > 0 there exists = (τ ) > 0, with (τ ) → 0 as τ → 0, such that for any initial condition a ∈ C n the solution of (3) satisfies In other words, for any control u(t) that is strictly positive for all t ≥ 0 we have that after any time τ > 0 all the normalized densities are strictly separated away from zero and from one. In particular, if the densities converge to a steady-state e then e i ∈ (0, 1) for all i. Contraction theory is a powerful tool for analyzing nonlinear dynamical systems [57, 58] , and has found applications in bio-molecular systems, control theory, synchronization of coupled nonlinear systems [59] , reaction-diffusion differential equations [60] , and more. For x ∈ R n , let |x| 1 := |x 1 | + · · · + |x n | denote the L 1 norm of x. For a non-singular matrix P ∈ R n×n , let |x| P,1 := |P x| 1 , i.e. the scaled L 1 norm of x. Proposition 4.2 Consider the RFMLK with a control u such that u(t) ≥ s > 0 for all t ≥ 0, and fix τ > 0. There exist a non-singular matrix P = P (τ ) and η = η(τ ) > 0 such that for any a, b ∈ C n , we have In other words, the RFMLK is contracting with respect to the scaled norm | · | P,1 after (the arbitrarily small) time delay τ . Proposition 4.2 implies several useful asymptotic properties of the RFMLK. These are described in the following subsections. In other words, the solution converges to e s for any initial condition, so the inital condition is "forgotten". The equilibrium e s represents a dynamic steady-state where the input and output flows from each site in the RFMLK are equal, and thus the densities in each site are constant. Angeli and Sontag [31] extended the notion of a monotone system to control systems. The next result shows that the RFMLK is a monotone control system. For two vectors v, w ∈ R n , we write v ≤ w if v i ≤ w i for all i = 1, . . . , n, and v w if v i < w i for all i = 1, . . . , n. Proposition 4.4 Fix two initial conditions a, b ∈ C n , with a ≤ b, and two controls u, v : Then the corresponding solutions of the RFMLK satisfy and In other words, if we consider two identical RFMLKs-the first with initial densities a i and the second with initial densities b i , with a i ≤ b i for all i, and apply a control u in the first and v in the second, with u(t) ≤ v(t) for all t ≥ 0-then at each time t ≥ 0 each density in the first RFMLK will be smaller or equal to the corresponding density in the second RFMLK. The next proposition analyzes the relation between the steady-state densities corresponding to constant control values. Then e s 1 ≤ e s 2 . We now turn to analyze the RFMLKN. For a matrix P ∈ R 1 × 2 , let P ∈ R 2 × 1 denote the transpose of P . Recall that every x i j ∈ [0, 1], and that the pool density satisfies z ∈ [0, ∞), so the state-space of the RFMLKN is For a ∈ Ω, let x(t, a) z(t, a) denote the solution of the RFMLKN at time t with the initial condition a. The next result follows immediately from the equations of the RFMLKN. The state space Ω in (9) is an invariant set for the dynamics of the RFMLKN that is, if a ∈ Ω then x(t, a) z(t, a) ∈ Ω for all t ≥ 0. In other words, every trajectory emanating from an initial condition in the state space remains in it for all t ≥ 0. The next result shows that trajectories that begin from an initial condition in Ω become uniformly separated from the boundary of Ω. . . , λ n + λ n−1 x n−1 + α n + β n u, For any x ∈ [0, 1] n all the entries of M (x) are nonnegative, so (5) is a cooperative dynamical system [54] . The matrix M (x) may become reducible for values x on the boundary of [0, 1] n . However, M (x) is irreducible for all x ∈ (0, 1) n . Thus, Proposition 4.7 guarantees that after an arbitrarily short time the matrix M (x(t)) and, thus J(x(t), u(t)), becomes an irreducible matrix. This will be used in analyzing the asymptotic properties of the RFMLKN described below. The next result shows that every level set contains a unique steady-state ribosome distribution in each mRNA and in the pool. The proof is based on the theory of monotone dynamical systems that admit a first integral, see [55, 56] . In other words, the RFMLKN admits a continuum of equilibrium points and any two solutions starting from initial conditions in the same level set of the system converge to the same equilibrium point. Thus, the rates λ i j , β i j , α i j and the total number of ribosomes s in the network determine a unique steady-state density profile in the RFMLKs and the pool. Equation (10) implies that for two initial conditions, the first one with a smaller total number of ribosomes than the second one, the corresponding equilibrium points e 1 and e 2 will be completely ordered: every steady-state density in e 1 will be strictly smaller than the corresponding density in e 2 . Example 7 To model a gene that is highly expressed with respect to other genes, consider an RFMLKN with a single RFMLK and a pool. Assume that the output function of the pool is G 1 (z) = z, and that the dimension of the RFMLK is n 1 = 2. The equations of the RFMLKN are thenẋ Any equilibrium point e = e 1 1 e 1 2 e z ∈ L s satisfies and e 1 1 + e 1 2 + e z = s. Various intracellular mechanisms may affect the parameters of the translation machinery. For example, the elongation rates depend on the interaction between the nascent peptide and the exit tunnel of the system [49] . Stress conditions increase ribosome abortion and drop-off. Due to competition for the finite pool of ribosomes, any change in the translation speed along a specific mRNA molecule will also indirectly affect the translation of other mRNAs in the network. The variability in the factors that affect translation in the cell requires models that can be used to Figure 9 : Trajectories of the RFMLKN in Example 7 for three initial conditions in L 1 . The unique equilibrium in L 1 is marked by a circle. analyze the sensitivity to parameter values. Another motivation for studying these issues comes from synthetic biology, for example, the recent interest in co-expression of multiple genes at a given, desired ratio [50, 51] . Our first result in this subsection analyzes the affect of a modification in the drop-off rate at one site of an mRNA molecule on the entire RFMLKN. We assume, without loss of generality, that the modification is in one of the rates in the first RFMLK. In other words, an increase in the drop-off rate in one of the RFMLKs yields an increase in the steady-state pool density and consequently an increase in the density in each site in all the other RFMLKs. The effect of increasing α 1 k on the steady-state in the first RFMLK is non-trivial. It is natural to expect a decrease in the density in each site of the first RFMLK. But, as more ribosomes accumulate in the pool the effective attachment rates in sites along the first RFMLK may also increase, leading to an increase in the density in certain sites. In general, the total effect on the first RFMLK will depend on all the parameters in the RFMLKN, and is thus difficult to predict. Example 8 Consider an RFMLKN with m = 2 RFMLKs of dimensions n 1 = 6 and n 2 = 3, parameters λ i 0 = 1, λ i j = 1, for all i, j, α 1 j = 0.01, for j = 1, 2, 4, 5, α 1 6 = 0, α 2 j = 0.01, for j = 1, 2, α 2 3 = 0, β 1 j = 0, for all j, β 2 1 = 0, β 2 j = 0.01, for j = 2, 3, and G i (z) = z, i = 1, 2. The initial condition is x i j = 0, for all i, j, and z(0) = 5. We simulated this RFMLKN for several values of the drop-off rate α 1 3 from the third site in RFMLK #1. Figures 10a and 10b show that increasing α 1 3 increases e 1 1 , and decreases e 1 i for all i > 1. As predicted in Thm. 4.2, it also increases e z , so more ribosomes accumulate in the pool leading to an increase in e 2 j for all j. Example 9 Consider the RFMLKN in Example 8, but now with β 1 4 = 2. Figure 11 shows that in this case an increase in the drop-off rate α 1 3 yields an increase in the steady-state values in the sites of RFMLK #1 located after the third site. This is because of an increase in the density of free ribosomes in the pool leading to more ribosomes attaching to the first RFMLK. The next result analyzes the "dual" case i.e. the affect of modifying one of the attachment rates in an RFMLK in the network. In other words, an increase in the attachment rate in one of the RFMLKs decreases the steadystate pool density and consequently decreases the density in each site in all the other RFMLKs. Example 10 Consider an RFMLKN with m = 2 RFMLKs of dimensions n 1 = 9 and n 2 = 3, and parameters λ i j = 1, for all i, j, α 1 j = 0.01, for all j except for α 1 9 = 0, α 2 j = 0.01, for j = 1, 2, α 2 3 = 0, β 1 j = 0 for all j except for β 1 3 , β 2 1 = 0, β 2 j = 0.01 for j = 2, 3, and with G i (z) = z, i = 1, 2. The initial condition is x i j = 0, for all i, j, and z(0) = 5. We simulated this RFMLKN for several values of the attachment rate β 1 3 . Figure 12b shows that as β 1 3 increases there is a decrease in the steady-state pool density (as more ribosomes bind to the third site of RFMLK #1) and consequently a decrease in the steady-state values in all the sites in RFMLK #2. As shown in Figure 12a , the effect of increasing β 1 3 on RFMLK #1 is non-trivial. The decrease in the pool density, decreases e 1 1 . However, the increase in the attachment rate in the third site leads to more ribosomes attaching to this site and consequently a higher density in sites 3, . . . , 9. Also, this creates a "traffic jam" along these sites and thus increases the density in site 2 as well. Our last result in this subsection analyzes the effect of modifying a hopping rate in one of the RFMLKs in the network. Consider an RFMLKN with m = 2 RFMLKs of dimensions n 1 = 9 and n 2 = 3, parameters λ i 0 = 2, λ i j = 2, for all i, j except for λ 1 5 , α 1 j = 0.1, for j ∈ {1, 2, . . . , 8}, α 1 9 = 0, α 2 j = 0.1, for j = 1, 2, α 2 3 = 0, β 1 j = 0.1 for all j except for β 1 1 = 0, β 2 1 = 0, β 2 j = 0.1 for j = 2, 3 and G i (z) = z, i = 1, 2. The initial condition is x i j = 0, for all i, j, and z(0) = 3. We simulated this RFMLKN for various values of the elongation rate λ 1 5 . Note that when λ 1 5 2 it is a bottleneck rate in RFMLK #1. Figure 13b shows that increasing λ 1 5 increases the pool density and thus the densities in all the sites along RFMLK #2. Figure 13a shows that increasing λ 1 5 yields a decrease in sites 3, 4, 5 in RFMLK #1, but an increase in all the other sites. Example 12 Consider again the RFMLKN in Example 11, but now the initial condition is x i j = 0, for all i, j, and z(0) = 10. In this case, Figure 14b shows that increasing the elongation rate λ 1 5 leads to a decrease in the steady-state pool density resulting in decreased steady-state densities in RFMLK #2. Figure 14a shows that in this case the effect of increasing λ 1 5 on RFMLK #1 is more intuitive: it decreases the densities in sites 1, . . . , 5 and increases the densities in sites 6, . . . , 9. We now analyze several additional mathematical properties of the RFMLKN. Recall that the dynamical systemẋ = f (x) is called cooperative if for any two initial condi- [54] . In other words, the flow preserves the (partial) ordering between the initial conditions. The next result shows that the RFMLKN is cooperative. If, furthermore, a = b then x(t, a) x(t, b) and z(t, a) < z(t, b), for all t ≥ 0. In other words, if we consider two initial conditions where the first corresponds to a smaller density in each site in each RFMLK and in the pool, then the corresponding solutions will satisfy the same relation at any time t ≥ 0. The next subsection shows that the flow of the RFMLKN is a non-expansive mapping. For a vector v ∈ R n , let |v| 1 = n i=1 |v i | denote the L 1 norm of v. In a contractive system, all solutions converge exponentially to one another. Since the RFMLKN admits more than one equilibrium, it is not a contractive system with respect to any norm. However, the next result shows that the L 1 distance between any two trajectories is non-expansive, i.e. it is bounded by the distance L 1 between the initial conditions. In particular, the difference between two "close" ribosomal density profiles will remain close for all t ≥ 0. In other words, the convergence to the equilibrium e Ls is monotone in the sense that the L 1 distance can only decrease with time. Many biological processes are excited by a periodic input. Proper functioning often requires entraining to the excitation, that is, converging to a periodic pattern with the same period as the excitation. A typical example is the ability of cells to coordinate their growth with the periodic cell-cycle division process. Translation seems to play an important role in this process. It is known for example that expression of the human translation initiation factor eIF3f peaks twice in the cell cycle: in the S and the M phases [64] . Entrainment is also important in the context of synthetic biology, for example, in designing a biological network that is coordinated by a single oscillator producing a common "clock signal" [68] . To study entrainment in the RFMLKN, assume that the parameters λ i j , α i j , β i j in all the RFM- LKs are not constant, but are time-varying functions, that are all jointly periodic with a period T > 0. More precisely, we assume that • There exists a (minimal) T > 0 such that all the non-negative time-varying rate functions λ i j (t), α i j (t) and β i j (t) are T-periodic. • There exists 0 < δ 1 ≤ δ 2 such that λ i j (t) ∈ [δ 1 , δ 2 ], for all i, j and all t ∈ [0, T ). We then refer to the network as the periodic RFM with Langmuir kinetics network (PRFMLKN). Note that a constant function is T -periodic for any T , so for example if one parameter in the network is T -periodic and all the others are constant then our assumptions hold. The next result shows that the PRFMLKN entrains. In other words, if the rates are T -periodic then all the densities in the mRNAs and the pool converge to a periodic pattern with period T , and thus so will the protein production rate in every mRNA. In particular, if a single parameter in one of the RFMLKs is T -periodic and all the other parameters are constant then the network entrains. Roughly speaking, this can be explained as follows. The T -periodic parameter will generate T -periodic variations in the pool density and this will generate T -periodic patterns in all the mRNA densities, as the pool feeds all the mRNAs. Again, this demonstrates the intricate coupling generated by the competition for shared resources. Example 13 Consider a network with m = 2 RFMLKs, of dimensions n 1 = 2 and n 2 = 3, and with G i (z) = tanh(z), i = 1, 2. All the rates in the network are equal to one, except for λ 2 2 (t) = 5 + 4 sin(2πt). Thus all the rates in the network are periodic with a common minimal period T = 1. The initial condition is z(0) = x i j (0) = 1/4 for all i, j. Figure 15 depicts the state-variables and the pool density as a function of t. Note that all the densities converge to a periodic pattern with period one. Note also that since the total number of ribosomes in conserved, maximal peaks in the density along the RFMLKs corresponds to minimal peaks in the pool density, and vice-versa. We derived and analyzed a novel and general network model of ribosome flow during largescale translation in the cell. This model encapsulates important cellular properties like ribosome drop-off, ribosome attachment at IRESs, and competition for a finite pool of free ribosomes. We analyzed the model using tools from systems and control theory, including contraction theory, and the theory of cooperative dynamical systems. The new model is an irreducible cooperative dynamical system admitting a first integral (the total density of ribosomes in the network). This implies that the system admits a continuum of linearly ordered equilibrium points, and that every trajectory converges to an equilibrium point. The system is also on the "verge of contraction" with respect to the L 1 norm. In addition, we proved that if one or more of the rates in the network are time-varying periodic functions with a common period T , then the densities along all the mRNAs and in the pool converge to a periodic solution with period T , i.e. the system entrains to a periodic excitation. An important question is the sensitivity of the network steady-state to variations in the mRNA parameters and the density of free ribosomes. We thoroughly analyzed this problem, and showed that a modification of a bio-physical property in a specific mRNA has two implications. what will be the effect on the densities and protein production rate in the mRNA that is modified, as this depends in a non-trivial way on the interactions between all the mRNAs and the pool. For example, an increase in the drop-off rate in a specific site in an mRNA may increase the pool density, thus increasing the attachment rates along this mRNA and leading to an increase in the density in some sites along this mRNA. These results highlight that analyzing the effect of any bio-physical property on translation in the cell must take into account the intricate effects of competition, especially when the competition for shared resources plays a major role, e.g. under stress conditions. We believe that the new model presented here provides a powerful framework for analyzing and re-engineering the translation process. One possible avenue for further research is in developing a quantitative and qualitative understanding of how viral mRNAs hijack the translation machinery and, in particular, whether the indirect effects of competition are enough to hamper the host's immune response. 6 Appendix: Proofs We begin by writing the RFMLK with input and output in the form: x where u : R + → R + is a scalar input function that takes non-negative values for any time t ≥ 0, y : R + → R + is a scalar non-negative output function, and Fix δ > 0. We will show that for any sufficiently small > 0 there exists K = K(δ, ) > 0 such that for every k ∈ {1, . . . , n} and every t ≥ 0 the condition For k = 1 the condition is simply x 1 ≤ , and theṅ Note that K 1 ≥ λ 0 s/2 > 0 for any > 0 sufficiently small. For k ∈ {2, . . . , n − 1} we havė Note that K k ≥ λ k−1 δ/2 > 0 for any > 0 sufficiently small. For k = n we havė Note that K n ≥ λ n−1 δ/2 > 0 for any > 0 sufficiently small. We conclude that (19) holds for K := min{K 1 , . . . , K n } ≥ min{λ 0 , . . . , λ n } min{s, δ}/2 > 0. By [36, Lemma 1] , this implies that for any τ > 0 there exists 1 = 1 (τ ) > 0, with 1 (τ ) → 0 as τ → 0, such that for any a ∈ C n the solution of (17) satisfies Let z i := 1 − x n+1−i , i = 1, . . . , n. Then (17) giveṡ It is not difficult to verify that this system also satisfies condition (19) , so by [36, Lemma 1] , for any τ > 0 there exists 2 = 2 (τ ) > 0, with 2 (τ ) → 0 as τ → 0, such that for any a ∈ C n the solution of (21) satisfies Let f denote the vector field in (17) , and let J := ∂ ∂x f denote its Jacobian. Then J = J 1 − diag(α 1 + β 1 u, , . . . , α n + β n u), where Note that J 1 is the Jacobian of an RFM with a time-varying initiation rate λ 0 u(t) ≥ λ 0 s, and that α i + β i u(t) ≥ 0 for all t. Note also that for any x ∈ int(C n ), all the entries in the superand sub-diagonal of J 1 are positive, so in particular J 1 (and thus J) is irreducible. Fix a, b ∈ C n and τ > 0. By Proposition 4.1, there exists = (τ ) > 0, such that It follows from (22) that J is a Metzler matrix, i.e. every off-diagonal entry of J is non-negative. Also, so every entry of K is non-negative. The results in [31] imply that the RFMLK is a monotone control system, so (7) holds. Now the definition of the output y implies (8) , and this completes the proof of Prop. 4.4. We already know that the limits e s 1 , e s 2 , and e := lim t→∞ x(t, a, v) exist. By monotonicity, x(t, a, u) ≤ x(t, a, v), for all t ≥ 0, and taking the limit as t → ∞ gives e s 1 ≤ e. Since the system is contractive, e = e s 2 , and this completes the proof of Prop. 4.5. For the sake of simplicity and to avoid cumbersome notation, we provide proofs of the theoretical results when m = 2, i.e. a network with two RFMLKs connected via a pool of free ribosomes (the proofs when m > 2 are very similar). We write the first RFMLK aṡ p n = λ n−1 p n−1 (1 − p n ) − λ n p n − α n p n + β n (1 − p n )u 1 , and the second asq The inputs to the RFMLKs are functions of the pool density and the pool dynamics iṡ The state-space of this network is Note that n i=1ṗ i + i=1q i +ż = 0, so is a first integral of the dynamics, that is, H(p(t), q(t), z(t)) ≡ H(p(0), q(0), z(0)). In other words, the total density of ribosomes in the network is conserved. For any s ≥ 0, we define the s level set of the first integral H by Thus, L s includes all the states in Ω with total density of ribosomes equal to s. Note that for s = 0, L 0 = {0} and the dynamics emanating from zero remains in zero for all time t ≥ 0. therefore, we will always consider L s with s > 0. We now restate and prove the persistence result in Prop. 4.7 for the case m = 2. Proposition 6.1 Consider the RFMLKN with m = 2. Fix s > 0. Then for any τ > 0 there exists (τ ) > 0, with (τ ) → 0 as τ → 0, such that for any initial condition in L s and any t ≥ τ the solution of the RFMLKN satisfies Proof: It is useful to denote p 0 := z, p n+1 = 0, and p −1 := p n . Using the fact that y 2 ≥ 0 We now show that the system with state-variables p 0 , . . . , p n satisfies the cyclic boundaryrepelling (CBR) condition in [33, Lemma 1] , that is, for any δ > 0 and any sufficiently small > 0, there exists K = K(δ, ) > 0 such that for each k = 0, . . . , n and each t ≥ 0 the condition p k (t) ≤ and p k−1 (t) ≥ δ implies thatẋ Indeed, for k = 0 we havė and since G i (0) = 0 and G i is continuous,ṗ 0 ≥ λ n δ/2 for all > 0 sufficiently small. For k ∈ {1, . . . , n} we havė It is also straightforward to verify that if p k (t) > 0 for some k ∈ {0, . . . , n} and t > 0 then p k (r) > 0 for all r ≥ t. It now follows from [33, Lemma 1] that for any τ > 0 there exists (τ ) > 0, with (τ ) → 0 as τ → 0, such that for any initial condition p 0 (0) . . . p n (0) Using a similar argument for the q system proves that for any t ≥ τ , Finally, arguing as in the proof of Prop. 4.1 completes the proof of Prop. 6.1. Our next goal is to prove the sensitivity results for the RFMLKN. It is useful to first write equations describing the steady-state of the RFMLK for a constant input u(t) ≡ v, with v > 0. By (3), Substituting the expressions for the f i s and g i s yields the following result. := λ n e n + n j=k+1 α j e j − n j=k+1 β j v(1 − e j ) λ k (1 − e k+1 ) , and v = w(e 1 , . . . , e n ) := λ n e n + n j=1 α j e j λ 0 (1 − e 1 ) + n j=1 β j (1 − e j ) . Note that the function w k is increasing in e k+1 , . . . , e n (and strictly increasing in e k+1 , e n ), and is decreasing in v. Also, the function w is increasing in e 1 , . . . , e n (and strictly increasing in e 1 , e n ). It is clear that for s = 0, L 0 = {e 0 }, with e 0 = 0, and all trajectories converge to e 0 . Fix s > 0. The Jacobian of the RFMLKN with m = 2 is where J pp is the Jacobian of an RFMLK with state-variables p 1 , . . . , p n , rates λ i , α i , β i and input u = G 1 (z) (see (22) ), J qq is the Jacobian of an RFMLK with state-variables q 1 , . . . , q , rates η i , γ i , δ i and input u = G 2 (z) (see (22) ), v := λ 0 (1 − p 1 )G 1 (z) + β 1 (1 − p 1 )G 1 (z) β 2 (1 − p 2 )G 1 (z) . . . β n (1 − p n )G 1 (z) T , w := η 0 (1 − q 1 )G 2 (z) + δ 1 (1 − q 1 )G 2 (z) δ 2 (1 − q 2 )G 2 (z) . . . δ (1 − q )G 2 (z) T , c := λ 0 G 1 (z) + α 1 + β 1 G 1 (z) α 2 + β 2 G 1 (z) . . . α n−1 + β n−1 G 1 (z) λ n + α n + β n G 1 (z) T , d := η 0 G 2 (z) + γ 1 + δ 1 G 2 (z) γ 2 + δ 2 G 2 (z) . . . γ −1 + δ −1 G 2 (z) η + γ + δ G 2 (z) T , and r : These equations imply that for any p q z T ∈ Ω the Jacobian matrix J is Metzler, so the RFMLKN is a cooperative dynamical system. Furthermore, for any p q z T ∈ int(Ω) all the entries on the super-and sub-diagonals of J pp , J qq are positive, and so are the first entry in v, w, and the first and last entry in c, d. This implies that for any p q z T ∈ int(Ω) the matrix J is irreducible. Combining this with Prop. 6.1 and the results in [56] on strongly cooperative dynamical systems with a first integral completes the proof of Theorem 4.1. We again prove for the special case m = 2, i.e. a network with two RFMLKs and a pool. Let e 1 i , i ∈ {1, . . . , n}, denote the steady-state in the first RFMLK, and e 2 j , j ∈ {1, . . . , }, denote the steady-state in the second RFMLK. Since the initial condition remains the same, we have: We prove that e z <ē z by contradiction. We consider two cases: e z =ē z ; and e z >ē z , and show that each of these cases yields a contradiction. Case 1. Assume that e z =ē z . In this case, there is no change in the input and parameter values in the second RFMLK, so e 2 j =ē 2 j for all j = 1, . . . , . Consider the first RFMLK. Since the steady-state input to this RFMLK remains the same, Prop. 7 in Ref. [35] , that states that increasing any of the detachment rates (without changing any other parameter) in the RFMLK decreases all the steady-state densities, implies thatē 1 j < e 1 j for all j. However, this contradicts (33) . In other words, after increasing α 1 k toᾱ 1 k , the steady-state input to each RFMLK is decreased. Thenē 2 j < e 2 j for all j = 1, . . . , . Combining this with (33) implies that e 1 j <ē 1 j for at least one index j. Let s ∈ {1, . . . , n} be the maximal index such that e 1 s <ē 1 s . Applying (30) inductively implies that e 1 j <ē 1 j for any j ∈ {s, s − 1, . . . , 1} (note that it follows from (30) that increasing α 1 k can only further increase the corresponding steady-state valueē 1 k ). Suppose that s < n. Then e 1 j <ē 1 j , j = 1, . . . , s, e 1 j ≥ē 1 j , j = s + 1, . . . , n. By (30), By the first equation in (29), G 1 (e z ) = λ 1 n e 1 n + n j=1 α 1 j e 1 j − n j=1 β 1 j G 1 (e z )(1 − e 1 j ) , and thus . Combining this with (34) and (37) implies that ifN ≥ N then G 1 (ē z ) ≥ G 1 (e z ), but this contradicts (34), so we conclude thatN < N . Now (39) implies thatD < D, i.e.ē 1 s+1 > e 1 s+1 . However, this contradicts (37), so we conclude that (36) cannot hold, that is, s = n. Applying (30) inductively implies that e 1 j <ē 1 j , j ∈ {n, n − 1, . . . , 1}. Eq. (31) gives G 1 (e z ) = λ 1 n e 1 n + n j=1 α 1 j e 1 , so G 1 (e z ) < G 1 (ē z ). This contradicts (34) , and we conclude that Case 2 is impossible. This completes the proof of Theorem 4.2. The proof is similar to the proof of Theorem 4.2 above and is thus omitted. In the first case, the input to each RFMLK satisfies G i (ē z ) = G i (e z ) Recall that the Jacobian J of the RFMLKN (with m = 2) is given in (32) Metzler, and a calculation shows that the sum of every column of J is zero. Hence, the L 1 matrix measure of J is zero Then f (t, z) = f (t+T, z) for all t ≥ 0 and z ∈ Ω, and H(x) Central dogma of molecular biology Molecular biology of the cell Cell and molecular biology The ribosome profiling strategy for monitoring translation in vivo by deep sequencing of ribosomeprotected mRNA fragments Genome-wide analysis of mRNA translation profiles in Saccharomyces cerevisiae The dynamics of supply and demand in mRNA translation Synthesis of proteins in Escherichia coli is limited by the concentration of free ribosomes: expression from reporter genes does not always reflect functional mRNA levels Quantitative assessment of ribosome drop-off in E. coli Novel mRNA-specific effects of ribosome drop-off on translation rate and polysome profile Mechanism of premature translation termination on a sense codon Tinkering with translation: protein synthesis in virus-infected cells Internal ribosome entry sites in eukaryotic mRNA molecules SARS-CoV-2 Nsp1 binds the ribosomal mRNA channel to inhibit translation Non-canonical translation in RNA viruses Exploiting internal ribosome entry sites in gene therapy vector design Mathematical and computational modelling of ribosomal movement and protein synthesis: an overview Kinetics of biopolymerization on nucleic acid templates Modeling translation in protein synthesis with TASEP: A tutorial and recent developments Totally asymmetric exclusion process with extended objects: a model for protein synthesis Nonequilibrium steady-states of matrix-product form: a solver's guide Totally asymmetric simple exclusion process with Langmuir kinetics Shock formation in an exclusion process with creation and annihilation Stochastic theory of protein synthesis and polysome: Ribosome profile on a single mRNA transcript Physics of transport and traffic phenomena in biology: from molecular motors and cells to organisms Nonequilibrium steady states of matrix-product form: a solver's guide Genome-scale analysis of translation elongation with a ribosome flow model Stability analysis of the ribosome flow model Ribosome Flow Model on a Ring Ribosome flow model with positive feedback Monotone control systems Ribosome flow model with different site sizes A model for competition for ribosomes in the cell Algorithms for ribosome traffic engineering and their potential in improving host cells' titer and growth rate A deterministic mathematical model for bidirectional excluded flow with Langmuir kinetics Entrainment to periodic initiation and transition rates in a computational model for gene translation Variability in mRNA translation: a random matrix theory approach Networks of ribosome flow models for modeling and analyzing intracellular traffic Strictly cooperative systems with a first integral Composite effects of gene determinants on the translation speed and density of ribosomes Global and local depletion of ternary complex limits translational elongation Far-from-equilibrium transport with constrained resources Translation with frameshifting of ribosome along mRNA transcript Preferential translation of internal ribosome entry sitecontaining mRNAs during the mitotic cycle in mammalian cells uORFs, reinitiation and alternative translation start sites in human mRNAs Internal ribosome entry sites in eukaryotic mRNA molecules Translational control in virus-infected cells Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems: an introduction to the theory of competitive and cooperative systems A comparative genomics study on the effect of individual amino acids on ribosome stalling Systematic comparison of 2A peptides for cloning multi-genes in a polycistronic vector Multicistronic IVT mRNA for simultaneous expression of multiple fluorescent proteins Production and application of multicistronic constructs for various human disease therapies Entrainment to periodic initiation and transition rates in a computational model for gene translation Monotone dynamical systems: an introduction to the theory of competitive and cooperative systems: an introduction to the theory of competitive and cooperative systems Strictly nonautonomous cooperative system with a first integral Cooperative irreducible systems of ordinary differential equations with first integral On contraction analysis for non-linear systems Contraction after small transients Global entrainment of transcriptional systems to periodic inputs Contraction methods for nonlinear systems: A brief introduction and some open problems Distributed nonlinear control algorithms for network consensus Oscillations in cell biology Translate to divide: control of the cell cycle by protein synthesis Expression of human eukaryotic initiation factor 3f oscillates with cell cycle in A549 cells and is essential for cell viability Comprehensive identification of cell cycle-regulated genes of the yeast Saccharomyces cerevisiae by microarray hybridization Genes adopt non-optimal codon usage to generate cell cycle-dependent oscillations in protein levels Formulation of the protein synthesis rate with sequence information Synthetic biology: applications come of age Models in systems biology: the parameter problem and the meanings of robustness