On the Ribosomal Density that Maximizes Protein Translation Rate

During mRNA translation, several ribosomes attach to the same mRNA molecule simultaneously translating it into a protein. This pipelining increases the protein translation rate. A natural and important question is what ribosomal density maximizes the protein translation rate. Using mathematical models of ribosome flow along both a linear and a circular mRNA molecules we prove that typically the steady-state protein translation rate is maximized when the ribosomal density is one half of the maximal possible density. We discuss the implications of our results to endogenous genes under natural cellular conditions and also to synthetic biology.


Introduction
The transformation of the genetic information in the DNA into functional proteins is called gene expression. Two important steps in gene expression are transcription of the DNA code into messenger RNA (mRNA) by RNA polymerase (RNAP), and then translation of the mRNA into proteins. During translation, complex macromolecules called ribosomes traverse the mRNA strand, decoding it codon by codon into a corresponding chain of amino-acids that is folded co-and post-translationally to become a functional protein [1]. The rate in which proteins are produced during the translation step is called the protein translation rate (or protein production rate).
According to current knowledge, translation takes place in all living organisms and under all conditions. Understanding the numerous factors that affect this dynamical process has important implications to many scientific disciplines including medicine, evolutionary biology, synthetic biology, and more.
Computational models of translation are becoming increasingly important as the amount of experimental findings related to translation rapidly increases (see, e.g. [2][3][4][5][6][7][8][9][10][11]). Such models are particularly important in the context of synthetic biology and biotechnology, as they can provide predictions on the qualitative and quantitative effects of various manipulations of the genetic machinery. Recent advances in measuring translation in real time [12][13][14][15] will probably further increase the interest in computational models that can integrate and explain the measured biological data. dynamics, as the total density is conserved. Both the RFM and RFMR are cooperative dynamical systems [32], but their dynamical properties turn out to be quite different [31].
The RFM [RFMR] has been applied to model and analyze ribosome flow along an open [circular] mRNA molecule during translation. Indeed, in eukaryotes the mRNA is often temporarily circularized (via non-covalent interactions), for example, by translation initiation factors [33]. Thus, a large fraction of the ribosomes that complete translating the mRNA reinitiate, and thus the RFMR is a good approximation of the translation dynamics in these circularized mRNAs. In addition, circular RNA forms (which include covalent RNA interactions) appear in all domains of life [34][35][36][37][38][39][40][41]. Specifically, it was recently suggested that circular RNAs can be translated in eukaryotes [39][40][41]. These cases are better approximated by the RFMR. Note that there are also cases of circular DNA [42,43]. This issue is not directly related to translation, yet transcription of such circular DNAs may be analyzed by a model similar to the RFMR [44] (see also [45]).
Here, we use the RFM [RFMR] to analyze the ribosomal density along a linear [circular] mRNA molecule that maximizes the steady-state protein translation rate. We refer to this density as the optimal density. This problem has already been studied before. For example, Zouridis and Hatzimanikatis [46] derived a deterministic, sequence-specific kinetic model for translation and studied the effect of the average ribosomal density on the steady-state translation rate. Their model assumes homogeneous elongation rates and open-boundary conditions, and includes all the elementary steps involved in the elongation cycle at every codon. Their simulations suggest that there exists a unique average density that corresponds to a maximal translation rate, see Figures 2A and 5A in [46] (see also [47]).
The RFM and RFMR are simpler models and thus allow to rigorously prove several analytic results on the optimal density. For a circular mRNA, we prove that there always exists a unique optimal density that maximizes the steady-state translation rate, and that it can be determined efficiently using a simple "hill climbing" algorithm. In addition, we show that under certain symmetry conditions on the rates the optimal density is one half of the maximal possible density.
In the case of a linear mRNA molecule, we prove that when the initiation and elongations rates are chosen to maximize the translation rate, under an affine constraint on the rates, the corresponding optimal density is one half of the maximal possible density (see Fig 2).
The remainder of this paper is organized as follows. The next section briefly reviews the RFM and the RFMR. Section 2 describes our main results. Section 3 summarizes the results, describes their biological implications, and suggests several directions for further research. The mathematical background and details are given in Appendix A. The proofs of all the results are placed in Appendix B.

The Ribosome Flow Model
Recall that the RFM models the ribosomes on the mRNA as "material" that flows between consecutive sites. For each site, the dynamics of the model can be expressed mathematically by the change in the amount of the material (i.e. ribosomes density) in that site, as a function of time, which is simply equal to the input rate to that site minus the output rate from that site. Since the model contains n sites, the RFM is expressed by n first-order ordinary differential  On the Ribosomal Density that Maximizes Protein Translation Rate equations describing the change in the amount of material in each site: . . .
where _ x i denote the change in the amount of material in site i as a function of time, i.e. _ x i ≔ d dt x i ðtÞ, i = 1, . . ., n. Note that this is a set of n nonlinear ordinary differential equations. If we define x 0 (t) ≔ 1 and x n+1 (t) ≔ 0 then Eq (1) can be written more succinctly as This can be explained as follows. The input rate to site i is . This flow is proportional to x i − 1 (t), i.e. it increases with the density at site i − 1, and to (1 − x i (t)), i.e. it decreases as site i becomes fuller. In particular, when this site is completely full, i.e. x i (t) = 1, there is no flow into this site. This corresponds to a "soft" version of a simple exclusion principle: the flow of particles into a site decreases as that site becomes fuller. Note that the maximal possible flow from site i − 1 to site i is the (i − 1)th transition rate λ i − 1 . Similarly, the output rate from site i, which is also the input rate to site i + 1, is given by The output rate of ribosomes from the chain is R(t) ≔ λ n x n (t), that is, the flow out of the last site. This is the number of proteins generated per time unit while translating the modeled mRNA template.
Suppose that we fix the values of the parameters in the RFM, that is, the transition rates. For any set of initial densities along the mRNA at the initial time zero, that is, the vector x(0) = [x 1 (0) . . . x n (0)] T , the RFM admits a unique dynamical solution x(t) = [x 1 (t) . . . x n (t)] T for any time t, where x i (t) is the density at site i at time t. Since every x i describes amount of normalized material in the sites, we always assume that their initial values are between zero and one. Combining this with the RFM dynamics, it can be shown that the amount of material in each site is bounded between zero and one for all time, i.e. x i (t) 2 [0,1], for all t ! 0 (see Section A.1). We denote the state-space, i.e. the space of all possible values of the state variables, by C n .
It was shown in [48] that the RFM is a cooperative dynamical system [32], and that this implies the following property. Every solution of the dynamics converges to a steady-state value, that is, the density x i (t) at site i converges, as time goes to infinity, to a specific value e i . We denote e ≔ [e 1 . . . e n ] T , i.e. the column vector of steady-state densities at each site. Furthermore, the value e will not depend on the initial value x(0). This means that if we simulate the RFM starting from any density of ribosomes on the mRNA the dynamics will always converge to the same final state (i.e., the same final ribosome density along the mRNA). Mathematically, the convergence property is written as lim t!1 x(t, x(0)) = e for all x(0) 2 C n (see also [49]). In particular, since x n (t) converges to e n , the translation rate R(t) = λ n x n (t) converges to the steady-state value: is analyzed below. The value of e does depend on the rates, i.e. e = e(λ 0 , . . ., λ n ) and satisfies e i 2 (0,1) for all i, that is, at steady-state each site is not completely empty nor completely full.

Ribosome Flow Model on a Ring
If we consider the RFM under the additional assumption that all the ribosomes leaving site n circulate back to site 1 then we obtain the RFMR. Just like the RFM, this is described by n nonlinear, first-order ordinary differential equations: . . .
Note that the only difference here relative to the RFM is in the equations describing the change of material in sites 1 and n. Specifically, the flow out of site n is the flow into site 1. This model assumes perfect re-cycling, but it also provides a good approximation when a large fraction of the ribosomes are re-cycled. Note that the model here is indifferent to the nature of the biochemical nature of the mRNA circularization; for example it can be both covalent or noncovalent RNA.
The RFMR can also be written succinctly as Eq (2), but now with every index interpreted modulo n. In particular, The total density of ribosomes along the chain at time t is given by i.e. the sum of the density at each site. Also, let s denote the total density of ribosomes along the chain at time 0, s ≔ H(x(0)). In the RFM, ribosomes enter and leave the chain and therefore H(x(t)) may vary with time t. In the RFMR, ribosomes that exit site n circulate back to site 1, so the total density is preserved for all time. This means that H(x(t)) = s for all t ! 0. The dynamics of the RFMR thus redistributes the particles between the sites, but without changing the total ribosome density. In the context of translation, this means that the total number of ribosomes on the (circular) mRNA is conserved. Mathematically, this means that H(x(t)) is a first integral of the RFMR (see Section A.2).
To understand the dynamics of the RFMR, we denote by L s the set of all possible ribosome density configurations such that the total density is equal to s. For example, L s includes the configuration where site 1 has a density s and all other sites are empty. It also includes the configuration where sites 1 and 2 have a density s/2 and all other sites are empty, etc. Mathematically, for s 2 [0,n], the s level set of H is where for p 2 R, p n denotes the column vector ½ p p . . . p T 2 R n .
Ref. [31] has shown that the RFMR is a strongly cooperative dynamical system, that every level set L s contains a unique equilibrium point e = e(s, λ 1 , . . ., λ n ), and that any trajectory of the RFMR emanating from any x(0) 2 L s converges to this equilibrium point. In particular, the translation rate converges to the steady-state value R = R(s, λ 1 , . . ., λ n ) (see Section A.2). This means the following. Fix an arbitrary value s 2 [0,n]. Then the RFMR, initiated with any configuration with total density s, will converge to the same final steady-state density e. This final density depends on the rates and the value s. For example if s = 0, corresponding to the initial condition x(0) = 0 n , then x(t) 0 n for all t ! 0, so e = 0 n . Similarly, s = n corresponds to the initial condition x(0) = 1 n and then clearly x(t) 1 n for all t ! 0, so e = 1 n . Since these two cases are trivial, below we will always assume that s 2 (0, n). In this case, e 2 Int(C n ), that is, the final density in each site will be strictly larger than zero and strictly smaller than one.
For more on the analysis of the RFM and the RFMR using tools from systems and control theory, see [31,[50][51][52][53][54]. For a general discussion on using systems and control theory in systems biology, see the excellent survey papers by Sontag [55,56].
The RFM models translation on a single isolated mRNA molecule. A network of RFMs, interconnected through a common pool of "free" ribosomes has been used to model simultaneous translation of several mRNA molecules while competing for the available ribosomes [57]. It is important to note that many analysis results for the RFM, RFMR, and networks of RFMs hold for any set of transition rates. This is in contrast to the analysis results on the TASEP model. Rigorous analysis of TASEP seems to be tractable only under the assumption that the internal hopping rates are all equal (i.e. the homogeneous case). In the context of translation, this models the case where all elongation rates are assumed to be equal.
The next section describes our main results on the optimal ribosome density.

Main Results
The average ribosomal density along the mRNA molecule at time t is simply the sum of all the site densities at time t, divided by the number of sites. We denote this average density by rðtÞ ≔ 1 n ðx 1 ðtÞ þ . . . þ x n ðtÞÞ. Recall that for every set of parameters in our models the state variables converge to a steady-state e. In particular, ρ(t) converges to the steady-state average ribosomal density: Note that since e i 2 [0,1] for all i, ρ 2 [0,1]. We are interested in analyzing the density that is obtained when the parameter values in the model are the ones that maximize the steady-state translation rate. Our results below show that there is a correspondence between the optimal translation rate and the optimal average ribosomal density: the optimal translation rate is obtained when and only when the average ribosomal density admits a specific value.

Optimal Density in the RFMR
Recall that in the RFMR the dynamical behavior depends on the transition rates and on the total density of ribosomes on the mRNA at time zero s ≔ P n i¼1 x i ð0Þ that, due to conservation of the total number of ribosomes, is equal to the total density at any time t. This means that the average ribosomal density is constant: Our first result shows that in the case of a circular mRNA there always exists a unique average density of ribosomes ρ Ã = s Ã /n that corresponds to a maximal steady-state translation rate of proteins. This means that in order to maximize the steady-state translation rate, the mRNA must be initialized with a total density s Ã (the distribution of this total density along the mRNA at time zero is not important). Initializing with either more or less than s Ã will decrease the steady-state translation rate with respect to the one obtained when the circular mRNA is initialized with total density s Ã . The optimal value s Ã depends on the transition rates. This is mathematically formulated as follows. Fix arbitrary transition rates λ i > 0, i = 1, . . ., n, and let R(s) ≔ R(s;λ 1 , . . ., λ n ) and e(s) ≔ e(s;λ 1 , . . ., λ n ) denote the steady-state translation rate and the steady-state ribosomal density profile, respectively, as a function of s.

Proposition 1
For any set of rates λ i > 0 in the RFMR there exists a unique value s Ã = s Ã (λ 1 , . . ., λ n ) 2 (0, n) that maximizes R(s). Furthermore, for this optimal value e Ã ≔ e(s Ã ) and and The optimality condition (5) can be explained as follows. If the total ribosome density s is very small then there will not be enough ribosomes on the circular mRNA and the translation rate will be small (for example, for s = 0 we have e = 0 n , and thus R = λ 1 e 1 (1 − e 2 ) = 0). In this case, the product of the e i s is small, so e 1 . . .e n < (1 − e 1 ). . .(1 − e n ) and Eq (5) does not hold. If s is very large traffic jams evolve on the mRNA and again the translation rate will be small (for example, for s = n we have e = 1 n , and thus R = λ 1 e 1 (1 − e 2 ) = 0). In this case, e 1 . . .e n > (1 − e 1 ). . .(1 − e n ) and Eq (5) does not hold. Thus, Eq (5) describes the unique point where the balance between too few and too many ribosomes on the circular mRNA molecule is optimal. The proof of this result (given in Appendix B) shows that the steady-state translation rate R(s) is strictly increasing on [0,s Ã ) and strictly decreasing on (s Ã , n], so a simple "hill climbing" algorithm can be used to find s Ã . The next example demonstrates Proposition 1 in the special case where all the elongation rates along the circular mRNA are identical. In this case, as shown in the example, the average ribosome density that maximizes the steady-state protein translation rate is exactly one half of the maximal possible total density. Example 1 Consider an RFMR with λ 1 = . . . = λ n , i.e. all the elongation rates are equal. Denote their common value by λ c . Then it follows from Eq (4) that 1 n c, c > 0, is an equilibrium point. In other words, any density profile with an identical density c at each site is a steady-state. Consider an initial condition with total density s. Then the trajectory satisfies x 1 (t) + . . . + x n (t) = s for all t ! 0. In particular, this must hold for the steady-state value, so nc = s or c = s/n. By uniqueness of the equilibrium point in every level set of H this implies that e = (s/n)1 n , so R = λ n e n (1 − e 1 ) = λ c (s/n)(1 − (s/n)). Thus, @R @s ¼ l c n 2 ðn À 2sÞ, so R(s) is strictly increasing [decreasing] on s 2 [0,n/2] [s 2 [n/2, n]] and therefore attains a unique maximum at s Ã = n/2. Then e Ã ≔ e(s Ã ) = (1/2)1 n and R Ã ≔ R(s Ã ) = λ c /4, and it is straightforward to verify that Eqs (5) and (6) hold. Note also that @ 2 R @s 2 ¼ À 2l c n 2 < 0, implying that R(s) is a strictly concave function.
The next example demonstrates the dependence of R(s) on s when the rates are not homogeneous.
Example 2 Consider an RFMR with dimension n = 3 and transition rates λ 1 = 2, λ 2 = 6, and λ 3 = 1/3. Note that e 3 is the maximal entry in e for all s. This is due to fact that the entry rate λ 2 = 6 into site 3 is high, and the exit rate λ 3 = 1/3 from site 3 is low.
In other words, in the typical case where the elongation rates along the circular mRNA are not all identical the average ribosome density that maximizes the protein translation rate is very close but not equal to half of the maximal possible total density. See Section A.3 for more explanations on Fig 4. For small values of n it is possible to provide more explicit results. The following result provides the maximal possible translation rate value, as a function of the elongation rates, in a very short RFMR with 2 or 3 sites.
Fact 1 For an RFMR with n = 2 the optimal values are s Ã = 1 and For an RFMR with n = 3 the optimal translation rate satisfies the equation: It is interesting to understand how the steady-state densities along the mRNA change as the total density along the circular mRNA varies slightly from the optimal value s Ã . This is important since in practice achieving a total number of ribosomes that equals exactly s Ã may be difficult, and the actual total number may be close, but not exactly equal to s Ã . For fixed transition rates, let e 0 i ≔ @ @s e i ðsÞ denote the change in the steady-state ribosome density at site i corresponding to a small change in the the total density s. We refer to e 0 i as the sensitivity of e i with respect to a change in the total density s.
The next result provides an expression for these sensitivities at the optimal density (i.e. the one that maximizes the protein translation rate). It can be used to compute how increasing/ decreasing the total number of ribosomes on the circular mRNA affects the optimal density profile.
Proposition 2 Consider an RFMR with dimension n. Fix rates λ i > 0, and let s Ã = s Ã (λ 1 , . . ., λ n ) and e Ã = e Ã (λ 1 , . . ., λ n ) be as defined in Proposition 1. Then ðe Ã Þ 0 ¼ v This means that if we change the density from s Ã to " s ≔ s Ã þ ε, where ε is very small, then the steady-state translation rate changes from R Ã to and substituting the numerical values yields Indeed, this agrees with the fact that the graph of R(s) attains a maximum at s Ã . On the Ribosomal Density that Maximizes Protein Translation Rate In Example 2 above the optimal value s Ã is close, but not equal to one half of the maximal possible total density n/2 = 3/2. The next result describes a specific case where the optimal total density is exactly one half of the maximal possible density.
Proposition 3 Suppose that the transition rates in the RFMR are symmetric, that is, Then s Ã = n/2 and e Ã i ¼ e Ã nþ1À i for all i. Thus, in this case the optimal average density is ρ Ã = (n/2)/n = 1/2, and the steady-state occupancies are also symmetric. In other words, when the elongation rates in the RFMR are symmetric the optimal average ribosome density is exactly one half of the maximal possible average ribosome density.
Note that condition (10) always holds for n = 2. Also, since a cyclic permutation of the rates leads to an RFMR with the same behavior, it is enough that Eq (10) holds for some cyclic permutation of the rates. For n = 3 this holds if at least two of the rates λ 1 , λ 2 , λ 3 are equal.
We note that a result similar to Proposition 3 is known for the homogeneous TASEP with periodic boundary conditions, i.e. that a loading of 50% maximizes the steady-state flow (see, for example, the fundamental diagram in Figure 4.1 in [28]).

Optimal Density in the RFM
Due to the open boundary conditions in the RFM, the number of particles along the chain as a function of time is not conserved. Thus, in this section we analyze the steady-state densities corresponding to the rates that yield a maximal steady-state translation rate. We refer to this density [rates] as the optimal density [rates].
Consider first the problem of maximizing the steady-state translation rate given a constraint on the weighted sum of the initiation and elongation rates. This is motivated by the fact that the biological resources are of course limited. For example, all tRNA molecules are transcribed by the same transcription factors (TFIIIB) and by RNA polymerase III. Hence, if the production of a specific tRNA is increased then the production of some other tRNA must decrease. This is captured by the constraint on the weighted sum of the rates, since any increase in one of the rates must be compensated by a decrease in some other rate. This was formulated in [51] by the following optimization problem.
Problem 1 Fix parameters b, w 0 , w 1 , . . ., w n > 0. Maximize R = R(λ 0 , . . ., λ n ), with respect to its parameters λ 0 , . . ., λ n , subject to the constraints: In other words, maximize the steady-state protein translation rate given an affine constraint on the initiation and elongation rates. Here b is the "total available biocellular budget" that can be spent on all the rates, and the positive values w i , i = 0, . . ., n, can be used to provide a different weighting to the different rates.
Problem 1 formalizes, using the RFM, an important problem in both synthetic biology and biotechnology, namely, determine the transition rates that maximize the protein translation rate, given the limited biomolecular budget. See Section A.4 for more details on Problem 1.
Here, our goal is to determine what is the steady-state density when the optimal rates are used (i.e. when we use the rates that maximize the steady-state translation rate under the constraints defined in Problem 1). We refer to this as the optimal density. Let e Ã i , i = 1, . . ., n, denote the steady-state density at site i corresponding to the optimal rates l Ã 0 ; . . . ; l Ã n in Problem 1.
The next example demonstrates that the optimal solution of Problem 1 typically corresponds to a ribosome density that is one half of the maximal possible ribosome density.
Example 4 Using a simple numerical algorithm we solved 10 5 instances of Problem 1 for an RFM with length n = 11 and total budget b = 1. In each instance the weights w i were drawn independently from a uniform distribution over the interval [0,1]. For each instance, we computed the optimal rates l Ã i and the corresponding average steady-state optimal density r Ã ≔ 1 n P n i¼1 e Ã i . Fig 5 depicts a normalized histogram (that is, the empirical probability) of the 10 5 values of ρ Ã . It may be observed that typically ρ Ã is close to 1/2. Similar results are obtained when the weights are drawn using other statistics, e.g. exponential, Rayleigh, and Gamma distributions.
In the case where all the weights in Problem 1 are equal we can also derive theoretical results on the structure of the optimal densities and the optimal average density.

Homogeneous Affine Constraint.
Recall that in Problem 1 the total weighted sum of the rates is bounded. The different weights can be used to provide a different "importance" for each rate. For example, if w 0 is much larger than the other weights then this means that any small increase in the initiation rate λ 0 will greatly increase the total weighted sum, thus typically forcing the optimal value l Ã 0 to be small. In this section we consider the case where all the weights are equal. This means that in the optimization problem all the rates have equal importance. We refer to this case as the homogeneous constraint case. Indeed, in this case the weights give equal preference to all the rates, so if the corresponding optimal solution satisfies l Ã i > l Ã j for some i, j then this implies that, in the context of maximizing R, λ i is "more important" than λ j . By Eq (20) (see Appendix A), we may assume in this case, without loss of generality, that w 0 = . . . = w n = b = 1, so the constraint in Problem 1 becomes The next result shows that under the homogeneous constraint (12) the steady-state densities corresponding to the optimal solution decrease along the chain, that is the steady-state density at site 1 is the largest, the density at site 2 is the second-largest, etc. It also implies that these steady-state densities are anti-symmetric with respect to the center of the chain. This property immediately implies that the optimal average density is one half (i.e. ρ Ã = 1/2). In other words, in this case the optimal solution to Problem 1 corresponds to a ribosome density that is one half of the maximal possible density.
Proposition 4 Consider Problem 1 with the homogeneous constraint (12). Then the optimal steady-state occupancies satisfy If n is even then and if n is odd then In both cases, the corresponding optimal density is ρ Ã = 1/2. The next example demonstrates the results in Proposition 4. Example 5 Consider Problem 1 for an RFM with n = 11 and the homogeneous constraint (12). Fig 6 depicts the optimal values l Ã i , i = 0, . . ., 11. It may be seen that the l Ã i s are symmetric, i.e. l Ã i ¼ l Ã 11À i , and that they increase towards the center of the chain. The corresponding steady-state distribution is e Ã = [0.5913, 0.5224, 0.5059, 0.5016, 0.5004, 0.5000, 0.4996, 0.4984, 0.4941, 0.4776, 0.4087] T (see Fig 7). It may be seen that the steady-state densities strictly decrease along the chain and are anti-symmetric with respect to the center of the chain.
Since the RFM is the dynamic mean-field approximation of TASEP, our results naturally raise the question of what is the optimal particle density in TASEP, that is, the density yielding the maximal possible output rate. Recall that each particle models a ribosome and each ribosome that leaves the chain produces a protein, so the output rate is also the protein translation rate. Rigorous analysis of this problem in TASEP seems to be non trivial.
We used a simple grid-search to address the problem of maximizing the steady-state flow in HTASEP (i.e. TASEP with all internal rates equal to one) with respect to the parameters α and β subject to the constraint w 1 α + w 2 β = b. For L = 11 and w 1 = w 2 = b = 1 the solution is α Ã = β Ã = 1/2, and the corresponding steady-state occupancies (computed using Eq (3.65) in [30]) are all equal to 1/2. Thus the average optimal density is ρ Ã = 1/2.
We also ran 10000 tests with w 1 and w 2 chosen from an independent uniform distribution on [0,1]. In each case, a simple grid-search was used to find the optimal rates. Fig 8 depicts a normalized histogram of the optimal average steady-state ribosome density in an HTASEP with L = 30. It may be seen that the typical optimal density is about 1/2. A similar result has been reported in [58] that used HTASEP with a superposition of open and periodic boundary conditions.
These simulation results corroborate the analytic results derived above for the RFM and RFMR.  (12). It may be seen that the optimal transition rates l Ã i s are symmetric around the center of the chain, with higher values at the center of the chain. doi:10.1371/journal.pone.0166481.g006

Discussion
A natural analogy for the cell is that of a factory operating complex and inter-dependent biosynthesis assembly processes [59]. Increasing the translation rate can be done by both operating several identical processes in parallel, and by pipelining every single process. In the context of translation, many mRNA copies of the same gene are translated in parallel, and the same transcript is simultaneously translated by several ribosomes. A natural question is what is the density of ribosomes along the transcript that leads to a maximal translation rate. It is clear that a very small density will not be optimal, and since the ribosomes interact and may jam each other, a very high density is also not optimal.
We studied this question using dynamical models for ribosome flow in both a linear and a circular mRNA molecule. Our results show that typically the optimal density is close to one half of the maximal possible density.
In synthetic biology and biotechnology optimizing the translation rate is a standard goal, and we believe that our results can provide guidelines for designing and re-engineering transcripts. However, in vivo biological regulation of mRNA translation may have several goals besides optimizing the translation rate. For endogenous genes there are many additional constraints that shape the transcript, translation rates, and ribosome densities. For example, it is known that evolution optimizes not only protein levels, but also attempts to minimize their production cost [60,61]. This cost may include for example the biocellular budget required for producing the ribosomes themselves. Thus, we do not expect that the protein levels of all genes will be maximal. Rather, we expect that translation is optimized for proteins that are required with high copy numbers (e.g. those related to house keeping genes and some structural genes).
Furthermore, it is important to mention that there are various additional constraints shaping the coding regions of endogenous genes. These include various regulatory signals related to various gene expression steps, co-translational folding, and the functionality of the protein [22,23,[62][63][64][65]. Thus, under these additional constraints we do not necessarily expect to see ribosome densities that maximize the translation rate. On the Ribosomal Density that Maximizes Protein Translation Rate Indeed, experimental studies of ribosome densities in various organisms demonstrate that on average 15% − 20% of the mRNA is occupied by ribosomes [66,67]. However, in 241 genes in S. cerevisiae more than 40% of the mRNA is occupied by ribosomes [66]. This suggests that a ribosome density that is close to 1/2 is frequent in certain specific mRNA molecules. In addition, it seems that under stress conditions ribosomal densities (and traffic jams) increase (see, e.g. [17]). Thus, under such conditions we expect more mRNAs with ribosome densities close to 1/2 (see, for example, [68]).
Interestingly, the reported results are also in agreement with genome-wide simulations of the RFM that were performed based on the modeling of all the endogenous genes of S. cerevisiae, as reported in [29]. Indeed, Fig 4C there shows the ribosome density, averaged over all the sites of all the mRNAs, as a function of the initiation rate. The maximal translation rate corresponds to an average density of about 1/2.
We note in passing that for an RFM with dimension n, with all the rates equal (i.e. λ 0 = . . . = λ n ), the average steady-state ribosomal density is 1/2 for all n, and that for an RFM with dimension n, λ 0 ! 1, and equal elongation rates (i.e. λ 1 = . . . = λ n ), the average steady-state ribosomal density is nþ1 2n , thus approaching 1/2 as n increases [69]. Further studies may consider optimizing the translation rate under various additional constraints. For example, it will be interesting to study the optimal ribosome density when taking into account also the biocellular cost of protein production, or under given constraints on the allowed density profile, etc. In addition, it will be interesting to study the optimal densities in more comprehensive models that include competition for the free ribosomes between several mRNA molecules [57]. Another important issue, that is not captured by the RFM and RFMR, is that every ribosome covers several codons. Developing and analyzing RFM/RFMR models with "extended objects" is an important challenge.
Finally, TASEP has been used to model and analyze many other natural and artificial processes including traffic flow and the movement of motor proteins. The problem of the optimal density is of importance in these applications as well.

A Appendix: Mathematical Background and Details
A.1 RFM Let x(t, a) denote the solution of Eq (1) at time t ! 0 for the initial condition x(0) = a. Since every state variable x i describes the occupation level at site i normalized between zero and one, this may also be interpreted as the probability that site i contains a ribosome. We always assume that a belongs to the closed n-dimensional unit cube: C n ≔ fx 2 R n : x i 2 ½0; 1; i ¼ 1; . . . ; ng: It is straightforward to verify that this implies that x(t, a) 2 C n for all t ! 0. In other words, C n is an invariant set of the dynamics [48]. This means that the occupancy levels always remain bounded between zero and one.
At steady-state, that is for x = e the left-hand side of all the equations in Eq (1) is zero, so l 0 ð1 À e 1 Þ ¼ l 1 e 1 ð1 À e 2 Þ ¼ l 2 e 2 ð1 À e 3 Þ This yields the following expressions for the final ribosome density profile: e n ¼ R=l n ; e nÀ 1 ¼ R=ðl nÀ 1 ð1 À e n ÞÞ; . . . e 2 ¼ R=ðl 2 ð1 À e 3 ÞÞ; and Combining Eqs (17) and (18) provides an elegant finite continued fraction [70] expression for the steady-state translation rate R in terms of the initiation and elongation rates along the coding region: Note that this equation admits several solutions for R, however, we are interested only in the unique feasible solution, i.e. the solution corresponding to e 2 Int(C n ). Indeed, this is the only solution that also yields a value between zero and one for the steady-state ribosome density at each site.
Note also that Eq (19) implies that Rðcl 0 ; . . . ; cl n Þ ¼ cRðl 0 ; . . . ; l n Þ; for all c > 0; ð20Þ that is, R(λ 0 , . . ., λ n ) is a homogeneous function of degree one. This means that if the initiation rate and all the translation elongation rates are multiplied by a factor c > 0 then the steadystate translation rate will also increase by a factor c.
Ref. [51] proved that R(λ 0 , . . ., λ n ) is a strictly concave function on R nþ1 þþ . Please see the figures in [51] for an intuitive explanation of this property.
HðxðtÞÞ ¼ Hðxð0ÞÞ; for all t ! 0: Let R = R(s, λ 1 , . . ., λ n ) denote the steady-state translation rate in the RFMR for any x(0)2L s . It follows from Eq (4) that R = λ i e i (1 − e i + 1 ), i = 1, . . ., n (recall that in the RFMR, every index is interpreted modulo n). It is straightforward to verify that for any c > 0 Rðs; cl 1 ; . . . ; cl n Þ ¼ cRðs; l 1 ; . . . ; l n Þ: ð22Þ This means that if the initiation rate and all the translation elongation rates are multiplied by a factor c > 0 then the protein translation rate will also increase by a factor c.
For any i 2 {1, . . ., n}, let k i ≔ λ 1 . . .λ i − 1 λ i + 1 . . .λ n , i.e. the product of all the elongation rates except for rate i, and let m ≔ P n i¼1 l i , that is, the sum of all the elongation rates. If s % 0 then all the e i s will be small, so we can ignore the terms 1 − e i in Eq (23), and this yields e i % k i s m , i = 1, . . ., n. A similar argument shows that if s % n then e i % 1 À k iÀ 1 ðnÀ sÞ m , i = 1, . . ., n. For the particular case in Example 2 this implies that when s % 0 then e % (s/22) [3 1 18] T . In particular, e 2 < e 1 < e 3 . When s % 3 then e % (s/22) [18s − 32 3s + 13 s + 19] T . In particular, e 1 < e 2 < e 3 . In other words, the relative ordering between the steady-state density at each site may change as the total density changes.
Regrading Proposition 2, it also implies the following. Indeed it follows from Eqs (9) and , i = 1, 2, . . ., n. This means that v i > v i − 1 if and only if 1 À e Ã i > e Ã iÀ 1 . In other words, the sensitivity to a change in the total optimal density at site i is larger than the sensitivity at site i − 1 if and only if the amount of "free space" at site i is larger than the occupancy at site i − 1.

A.4 Optimal Density in the RFM
Consider Problem 1. It has been shown in [51] that the optimal solution l Ã 0 ; . . . ; l Ã n always satisfies P n i¼0 w i l Ã i ¼ b. Of course, by scaling the w i s we may always assume that b = 1. Combining this with the strict concavity of the steady-state translation rate R(λ 0 , . . ., λ n ) in the RFM implies that Problem 1 is a convex optimization problem that admits a unique optimal solution l Ã 2 R nþ1 þþ . Also, this solution can be determined efficiently using numerical algorithms that scale well with n.

B Appendix: Proofs
Proof of Proposition 1. It follows from known results on the solutions of ODEs that e i is continuous in s for all i. It is known that every e i is strictly increasing in s (see [Theorem 1] in [31] Multiplying these n inequalities, and using the fact that e 0 To prove the converse implication, assume that Eq (25)  where a i ≔ e 0 i ð1 À e iþ1 Þ, and b i ≔ e i e 0 iþ1 . This means that a ℓ > b ℓ for some index ℓ 2 {1, . . ., n}. Since R 0 = λ ℓ (a ℓ − b ℓ ), it follows that R 0 > 0. Thus, we showed that R 0 > 0 if and only if Q n i¼1 ð1 À e i Þ > Q n i¼1 e i . The proof that R 0 < 0 if and only if Q n i¼1 ð1 À e i Þ < Q n i¼1 e i is similar. This implies that R 0 = 0 if and only if Q n i¼1 ð1 À e i Þ ¼ Q n i¼1 e i , and this completes the proof of Proposition 5.
We can now complete the proof of Proposition 1. Let pðsÞ ≔ Q n i¼1 ð1 À e i Þ, and qðsÞ ≔ Q n i¼1 e i . Then p(0) = 1, p(n) = 0, q(0) = 0, and q(n) = 1. The strict monotonicity of every e i implies that p(s) [q(s)] is a strictly decreasing [increasing] function in the interval s 2 [0,n]. This implies that there is a unique s Ã 2 [0,n] such that p(s Ã ) = q(s Ã ). By Proposition 5, this is the unique maximizer of R(s), and for s = s Ã : e Ã 1 . . . e Ã n ¼ ð1 À e Ã 1 Þ . . . ð1 À e Ã n Þ: Also, Space, and the US-Israel Binational Science Foundation. The research of MM is also supported by a research grant from the Israel Science Foundation.