Maximizing Protein Translation Rate in the Nonhomogeneous Ribosome Flow Model: A Convex Optimization Approach

Translation is an important stage in gene expression. During this stage, macro-molecules called ribosomes travel along the mRNA strand linking amino-acids together in a specific order to create a functioning protein. An important question is how to maximize protein production. Indeed, translation is known to consume most of the cell's energy and it is natural to assume that evolution shaped this process so that it maximizes the protein production rate. If this is indeed so then one can estimate various parameters of the translation machinery by solving an appropriate mathematical optimization problem. The same problem also arises in the context of synthetic biology, namely, re-engineer heterologous genes in order to maximize their translation rate in a host organism. We consider the problem of maximizing the protein production rate using a computational model for translation-elongation called the ribosome flow model (RFM). This model describes the flow of the ribosomes along an mRNA chain of length n using a set of n first-order nonlinear ODEs. It also includes n+1 positive parameters: the ribosomal initiation rate into the mRNA chain, and n elongation rates along the chain sites. We show that the steady-state translation rate in the RFM is a strictly concave function of its parameters. This means that the problem of maximizing the translation rate under a suitable constraint always admits a unique solution, and that this solution can be determined using highly-efficient algorithms for solving convex optimization problems even for large values of n. Furthermore, our analysis shows that the optimal translation rate can be computed based only on the optimal initiation rate and the elongation rate of the codons near the beginning of the ORF. We discuss some applications of the theoretical results to synthetic biology, molecular evolution, and functional genomics.

Abstract. Translation is an important stage in gene expression. During this stage, macro-molecules called ribosomes travel along the mRNA strand linking amino-acids together in a specific order to create a functioning protein.
An important question, related to many biomedical disciplines, is how to maximize protein production. Indeed, translation is known to consume most of the cell's energy and it is natural to assume that evolution shaped this process so that it maximizes the protein production rate. If this is indeed so then one can estimate various parameters of the translation machinery by solving an appropriate mathematical optimization problem. The same problem also arises in the context of synthetic biology, namely, re-engineer heterologous genes in order to maximize their translation rate in a host organism. We consider the problem of maximizing the protein production rate using a computational model for translation-elongation called the ribosome flow model (RFM). This model describes the flow of the ribosomes along an mRNA chain of length n using a set of n first-order nonlinear ordinary differential equations. It also includes n + 1 positive parameters: the ribosomal initiation rate into the mRNA chain, and n elongation rates along the chain sites. We show that the steady-state translation rate in the RFM is a strictly concave function of its parameters. This means that the problem of maximizing the translation rate under a suitable constraint always admits a unique solution, and that this solution can be determined using highly-efficient algorithms for solving convex optimization problems even for large values of n. Furthermore, our analysis shows that the optimal translation rate can be computed based only on the optimal initiation rate and the elongation rate of the codons near the beginning of the ORF. We discuss some applications of the theoretical results to synthetic biology, molecular evolution, and functional genomics.
to the 3 ′ end [1]. Each triplet of the mRNA consecutive nucleotides, called a codon, is decoded by a ribosome into a corresponding amino-acid. The rate in which proteins are produced during the translation step is referred to as the protein production rate or translation rate.
The translation process occurs in all organisms, in almost all cells, and in almost all conditions. Thus, understanding translation has important implications in many scientific disciplines, including medicine, biotechnology, functional genomics, evolutionary biology, and more. The amount of biological findings related to translation increases at an exponential rate and this leads to considerable interest in computational models that can integrate and analyze these findings (see, e.g., [68,12,24,36,62,61,10,56,15,49]).
A fundamental challenge in biotechnology and synthetic biology is to control the expression of heterologous genes in a host organism in order to synthesize new proteins or to improve certain aspects of the host fitness [53,43,5]. Computational models of translation are also important in this context, as they allow one to simulate and analyze the effect of various manipulations of the genomic machinery on the translation process.
A conventional computational model of translation-elongation is the totally asymmetric simple exclusion process (TASEP) [57,69]. TASEP is a stochastic model that describes particles moving along a one-dimensional lattice of sites. The term totally asymmetric is used to indicate unidirectional motion along the chain. Each site can be either empty or occupied by a single particle. This captures interaction between the particles, as a particle in site i blocks the movement of a particle in site i − 1. Hence, the term simple exclusion. At each time instant, the sites are scanned and provided that a site is occupied by a particle and the next site is empty, the particle hops to the next site with some probability. The two sides of the chain are connected to particle reservoirs, and particles can hop into the chain (if the first site is empty) and out of the chain (if the last site is full). TASEP is a fundamental model in nonequilibrium statistical mechanics that has been used to model numerous natural and artificial processes [55]. Analysis of TASEP is based on determining the probabilities of steady-state configurations using matrix products (see the excellent review paper [7]). The ribosome flow model (RFM) [51] is a deterministic model for translationelongation that can be obtained via a mean-field approximation of TASEP (see, e.g., [55, section 4.9.7] and [7, p. R345]). The RFM for a chain with n sites includes n first-order, nonlinear ordinary differential equations and n + 1 positive parameters: the initiation rate λ 0 , and elongation rates λ i , i = 1, 2, . . . , n, between every two consecutive sites.
There are indications that in some genes all the elongation rates along the mRNA chain are (approximately) equal [29]. This may be modeled by assuming constant elongation rates in the RFM. This yields the homogeneous ribosome flow model (HRFM) [39] that includes only two positive parameters: the initiation rate λ 0 and the constant elongation rate λ c .
In a previous study [67], we have shown that the steady-state protein translation rate in the HRFM, denoted R = R(λ 0 , λ c ), is a concave function of the parameters λ 0 , λ c . The proof of this result is based on analyzing the Hessian matrix H of R in the HRFM. Note that H has dimensions 2 × 2 for all n. However, the assumption of equal elongation rates is often too strong. For example, it was shown that factors such as the adaptation of codons to the tRNA pool [13,33,11], folding of the mRNA [62,11], and local amino acid charge [62,11,9] affect translation elongation speed. This induces variations between different elongation rates. In these cases, the HRFM is not a suitable model, and one must use the RFM. The steady-state translation rate in the RFM is a function of n + 1 parameters, i.e., R = R(λ 0 , . . . , λ n ). In this paper, we show that R(λ 0 , . . . , λ n ), is a strictly concave function of its n + 1 positive parameters. Here the Hessian matrix has dimensions (n + 1) × (n + 1), and it seems that the approach applied in [67] cannot be extended to handle the RFM. The proof of our main result is thus based on an entirely new technique.
To explain the importance of the strict concavity of R, consider Fig. 1 that depicts, for simplicity, a scalar strictly concave function y = f (x). Strict concavity in this case means the following. Given any two different values x 1 , x 2 , with corresponding function values y 1 = f (x 1 ) and y 2 = f (x 2 ), let l = l(x) denote the line that connects the points (x 1 , y 1 ) and (x 2 , y 2 ). Then f (x) > l(x), for all x ∈ (x 1 , x 2 ). In other words, the graph of the function lies above the line l(x).
Concave functions have many useful and desirable properties. First, a concave function is differentiable almost everywhere. Second, recall that a point x m is called a local maximum of a function f if the function values in some neighborhood of x m are smaller than or equal to f (x m ). It is a global maximum if the function values in its entire domain of definition are smaller than or equal to f (x m ). For a concave function, any local maximum is also a global maximum. If the function is strictly concave then this maximum is unique.
Furthermore, strict concavity implies that a simple "hill climbing" algorithm can be used to find the global maximum. In the depicted one dimensional function, this can be explained as follows. Select an arbitrary point x 0 in the domain of definition of f as a candidate for a maximum point. Next, determine two points x − 0 and x + 0 that are "close" to x 0 and satisfy . If y 0 ≥ y − 0 and y 0 ≥ y + 0 then x 0 is a local, and thus global, maximum of the function and the algorithm terminates. Otherwise, at least one of the two values y − 0 , y + 0 is larger than y 0 . The corresponding point, i.e., x − 0 or x + 0 , becomes the new candidate for a maximum, and the algorithm is iterated. Under mild assumptions, this simple algorithm is guaranteed to converge to the global maximum of the concave function. More generally, there exist highly-efficient algorithms for finding the global maximum of multi-dimensional concave functions [8].
A function g is called [strictly] convex if −g is [strictly] concave. Thus, the problem of finding the maximum value of a concave function is equivalent to the problem of finding the minimum value of a convex function. A famous quote by R. T. Rockafellar states that: "...the great watershed in optimization isn't between linearity and nonlinearity, but convexity and nonconvexity." [52] We note in passing that a linear function y(x) = ax + b is both concave and convex.
Summarizing, our main result implies that the problem of maximizing the protein translation rate, under a simple constraint on the RFM parameter values, admits a unique solution, and that this solution can be found numerically using highlyefficient algorithms. It is important to note that many systems and processes have been modeled and analyzed using TASEP. These include translation, traffic flow, molecular motors, surface growth, the movement of ants, and more [55]. All these processes may also be modeled using the RFM, and the problem of maximizing R seems to be of importance in all of them.
We now describe some possible applications of the main result in the context of translation. A recent work [20] studied the effect of the intracellular translation factor abundance on the protein production rate. Abundance of the encoded translation factor was experimentally manipulated to a sub-wild-type level [20] using the tet07 construct. The reported results suggest that the mapping from levels of translation factors to protein production rate is concave (see Fig. 1 in [20]). This may provide an experimental support to the results presented in this paper. Note that [20] used the model organism S. cerevisiae that is known to have non-constant elongation rates [62,13]. Thus, the RFM, and not the HRFM [67], is a better computational model for describing these experiments.
Translation is known to consume most of the cell's energy [47,60,1]. A reasonable assumption is that in organisms under strong evolutionary pressure the genomic machinery has evolved so that it optimizes the translation rate given the available resources. This assumption can be studied in the context of the RFM since the concavity of the translation rate implies that one can easily determine the optimal parameter values, and then compare them to biological findings. This may help in understanding the level of selection pressure acting on the genomes of various organisms and the evolutionary changes in various micro-organisms [16].
In synthetic biology, an important problem is to re-engineer a genetic system by manipulating the transcript sequence, and possibly other intra-cellular variables, in order to obtain an optimal translation rate. Using our results on the RFM can provide verifiable predictions on how this can be done efficiently. Another related problem is optimizing the translation efficiency and protein levels of heterologous genes in a new host [47,60,23,30]. These genes actually compete with endogenous genes for the available resources, e.g., initiation factors. Consuming too much resources by the heterologous gene may kill the host [47,60]. Thus, any optimization of the protein translation rate should not consume too many resources, as otherwise the fitness of the host may be significantly reduced. This seems to fit well with the constrained optimization problem that we pose here for the RFM.
The remainder of this paper is organized as follows. Section 2 briefly reviews the RFM. Section 3 presents the main results. Section 4 describes the implications of our results to systems biology, evolution, and synthetic biology, and describes several possible directions for further research. To streamline the presentation, all the proofs are placed in the Appendix.

Preliminaries
The RFM [51] is a deterministic mathematical model for translation-elongation. In the RFM, mRNA molecules are coarse-grained into a unidirectional chain of n sites of codons. The RFM is a set of n first-order nonlinear ordinary differential equations: . . . Here, implies that site i is completely empty [completely full] at time t. The parameter λ 0 > 0 is the initiation rate into the chain, and λ i > 0, i ∈ {1, .., n}, is a parameter that controls the flow from site i to site i + 1. In particular, λ n controls the output rate at the end of the chain. 5 The rate of ribosome flow into the system is λ 0 (1 − x 1 (t)). The rate of ribosome flow exiting the last site, i.e., the protein production rate, is λ n x n (t). The rate of ribosome flow from site i to site i + 1 is Fig. 2). Note that this rate increases with x i (t) (i.e., when site i is fuller) and decreases with x i+1 (t) (i.e., when the consecutive site is becoming fuller). In this way, the RFM, just like TASEP, takes into account the interaction between the ribosomes in consecutive sites.
We emphasize that in the RFM the state-variables take values in the closed interval [0, 1] and are not limited to the values {0, 1}. This is different from TASEP, where a site can be either empty or full. Indeed, the x i s in the RFM may be interpreted as time-averaged occupancy levels in TASEP, and this average takes values in [0, 1]. Let x(t, a) denote the solution of (1) at time t ≥ 0 for the initial condition x(0) = a. Since the state-variables correspond to normalized occupation levels, we always assume that a belongs to the closed n-dimensional unit cube: It is straightforward to verify that this implies that x(t, a) ∈ C n for all t ≥ 0. In other words, C n is an invariant set of the dynamics [40].
Let int(C n ) denote the interior of C n . It was shown in [40] that the RFM is a monotone dynamical system [58] and that this implies that (1) admits a unique equilibrium point e ∈ int(C n ). Furthermore, lim t→∞ x(t, a) = e, for all a ∈ C n .
This means that all trajectories converge to the steady-state e.
We note in passing that monotone dynamical systems have recently found many applications in systems biology, see, e.g., [3,34,59] and the references therein.
For x = e, the left-hand side of all the equations in (1) is zero, so Denoting the steady-state translation rate by yields where we define e n+1 := 0. Also, . . . and Combining (5) and (6) provides a finite continued fraction [35] expression for R: Note that this equation has several solutions for R (and thus also several solutions for e n = R/λ n ), however, we are interested only in the unique feasible solution, i.e. the solution corresponding to e ∈ int(C n ).
Eq. (7) may be written as p(R) = 0, where p is a polynomial of degree ⌈(n + 1)/2⌉ in R with coefficients that are algebraic functions of the λ i s. For example, for n = 3, (7) yields Recent biological findings suggest that in some cases the transition rate along the mRNA chain is approximately constant [29]. This may be also the case for gene transcription [17]. To model this, Ref. [39] has considered the RFM in the special case where that is, the transition rates λ i , i = 1, 2, . . . , n, are all equal, and λ c denotes their common value. Since this homogeneous ribosome flow model (HRFM) includes only two parameters, λ 0 and λ c , the analysis is simplified. In particular, (7) becomes where λ c appears a total of n times.
Several recent papers analyzed the RFM or HRFM. To model ribosome recycling (see, e.g., [44] and the references therein), Ref. [41] has considered a closed-loop RFM with a positive linear feedback from the output R to the input λ 0 . It has been shown that the closed-loop system admits a unique globally asymptotically stable equilibrium point. In [38], it has been shown that the state-variables (and thus the protein production rate) in the RFM entrain to periodically time-varying initiation and/or transition rates. This provides a computational framework for studying entrainment to a periodic excitation (e.g., the cell cycle) at the translation level. The HRFM with an infinitely-long chain, (i.e. with n → ∞) was considered in Ref. [66]. There, a simple closed-form expression for e ∞ := lim n→∞ e n was derived, as well as explicit bounds for |e ∞ − e n | for all n ≥ 2.
In the RFM the steady-state production rate R is a function of the positive parameters λ 0 , . . . , λ n . In this paper, we study the dependence of R on these parameters. Our results are based on a novel, linear-algebraic approach linking the protein translation rate to the maximum eigenvalue of a symmetric, non-negative tridiagonal matrix whose components are functions of the λ i s.

Concavity
The next result is the main result in this section. Recall that all the proofs are placed in the Appendix. Let R n+1 Theorem 1 Consider the RFM with dimension n. The steady-state translation rate R = R(λ 0 , . . . , λ n ) is a strictly concave function on int(R n+1 + ). The next example demonstrates Theorem 1.
Example 1 Consider the RFM with n = 1. In this case, (5) and (6) yield e 1 = R/λ 1 and e 1 = 1 − R/λ 0 , so Fig. 3 depicts R(λ 0 , λ 1 ) as a function of its arguments. It may be seen that this is indeed a strictly concave function on int(R 2 + ). Recall that a function f : R k + → R is called positively homogeneous of degree m if f (cx) = c m f (x) for all c > 0 and all x ∈ R k + . For example, the function f (x 1 , x 2 ) := x 1 x 2 + x 2 2 is positively homogeneous of degree 2. The following result follows immediately from the fact that R always appears in (7) only in terms of the form R/λ i .  In other words, R(cλ 0 , . . . , cλ n ) = cR(λ 0 , . . . , λ n ), for all c > 0.
From a biophysical point of view this means that multiplying the initiation rate and all the elongation rates by the same factor c > 0 increases of the steady-state production rate by a factor of c. This also means that the steady-state occupancy levels e i , i = 1, 2, . . . , n, remain unchanged with respect to such a multiplication.

Constrained Maximization of the Protein Translation Rate
Consider the problem of determining the parameter values λ 0 , . . . , λ n that maximize R (or, equivalently, that minimize −R) in the RFM. Obviously, to make this problem meaningful we must constrain the possible parameter values. This leads to the following constrained optimization problem.
In other words, given an affine constraint on the total rates, namely, the initiation rate λ 0 and the transition rates λ 1 , . . . , λ n , maximize the protein translation rate. The constraint on λ i , i = 0, 1, . . . , n, may be related to factors such as the abundance of intracellular ribosomes, initiation factors, intracellular tRNA molecules and elongation factors. The values w i , i = 0, 1, . . . , n, can be used to provide different weighting to the different rates. It is not difficult to show that the optimal solution λ * of Problem 1 always satisfies λ * ∈ int(R n+1 + ). Theorem 1 implies that Problem 1 is a convex optimization problem [8]. It thus enjoys many desirable properties.
The next result shows that increasing any of the λ i s increases the translation rate.
In other words, increasing either the initiation rate or the elongation rate at any site improves the production rate.
Remark 1 Proposition 1 implies that the first constraint in (12) can be replaced by n i=0 w i λ i = b.
It is clear from (11) that in general an algebraic expression for R in terms of λ 0 , . . . , λ n does not exist. It is possible however to give an algebraic expression for the maximal value R * as a function of just two optimal parameter values, namely, λ * 0 and λ * 1 , and the parameters in the affine constraint.

Theorem 2 Consider Problem 1 for the RFM with dimension n. Then
In other words, the optimal translation rate R * can be computed given the optimal initiation rate and the first optimal elongation rate (and their corresponding weights in the affine constraint). This result holds regardless of the length of the transcript.
It is interesting to note that several biological studies showed that various signals encoded in the 5'UTR and the beginning of the ORF can predict the protein levels of endogenous genes with relatively high accuracy [33,32,70,63,47]. Maximization with equal constraint weights It is interesting to consider the specific case where all the weights w i in the constrained optimization problem are equal. Indeed, in this case the weights give equal preference to all the rates, so if the optimal solution satisfies λ * i > λ * j for some i, j then this may be interpreted as saying that, in the context of maximizing R, λ i is "more important" than λ j . Fig. 5 depicts the optimal values λ i for the case where b = 1 and w i = 1 for all i. In other words, the constraint is n i=0 λ i = 1. Three cases are shown corresponding to n = 30, n = 10, and n = 4. The optimal values were found numerically using a simple search algorithm that is guaranteed to converge for convex optimization problems.
It may be observed that the optimal transition rates are symmetric with respect to the index i = n/2. In general, the transition rate λ * n/2 is larger than all other rates and the optimal values decrease as we move towards any edge of the chain. The difference between λ * n/2 and λ * 0 (or λ * n ) is always visible, but the difference between λ * n/2 and λ * (n/2)±i , with i small, becomes negligible as n increases. Intuitively, these results may be interpreted as follows. The importance of an elongation rate (or the corresponding site) depends on its "centrality", or the mean distance of this site to other sites in the chain. Site n/2 is thus always the most "important" site in the chain. As n increases, the sites near the middle site have almost the same mean distance to the other sites, and thus become almost as important. Fig. 6 depicts the optimal translation rate R * as a function of n for two different constraints: n i=0 λ i = n and n i=0 λ i = n 1.03 . The first case corresponds to the scenario where the total available resources increases linearly with n (i.e., b = n). It may be observed that in this case the optimal translation rate R * decreases monotonically with n. On the other hand, increasing the total available resources by a rate which is slightly larger than a linear rate (i.e., b = n 1.03 ) changes the behavior; R * in this case increases monotonically with n. This result suggests that in order to maintain the same optimal translation rate value as n increases, the total allocated resources should increase at a rate that is slightly higher than a linear rate in n.

Discussion
The RFM is a deterministic mathematical model for translation-elongation. It can be derived via a mean-field approximation of a fundamental model from non-equilibrium statistical mechanics called TASEP. The RFM encapsulates both the simple exclusion and the total asymmetry properties of the stochastic TASEP model. The RFM is characterized by an order n, corresponding to the number of sites along the mRNA strand, a positive initiation rate λ 0 and a set of positive alongation rates λ 1 , . . . , λ n .
In this paper, we show that the steady-state protein translation rate R = R(λ 0 , . . . , λ n ) in the RFM is a strictly concave function of its (positive) parameters. This implies that: (1) a local maximum of R is the global maximum (and this maximum is unique); and (2) the problem of maximizing the steady-state protein translation rate under an affine constraint on the RFM parameters is a convex optimization problem. Such problems can be solved numerically using highly-efficient algorithms. The constraint here aims to capture the limited biosynthetic budget of the cell. We now describe the possible implications of these results in various disciplines including biology, synthetic biology, molecular evolution, and functional genomics. As mentioned above, the functional dependence of the translation rate on various variables can also be examined experimentally. A recent paper [20] studied the effect of the intracellular translation factor abundance on protein synthesis. Experiments based on a tet07 construct were used to manipulate the production of the encoded translation factor to a sub-wild-type level, and measure the translation rate (or protein levels) for each level of the translation factor(s). An analysis of Fig. 1 in [20] suggests that the mapping from levels of translation factors to the translation rate is indeed concave. Our results thus provide the first mathematical support for the observed concavity in these experiments.
In synthetic biology, re-engineering gene expression is frequently used to synthesize proteins for medical and agricultural goals [53,43,5,21,42]. For example, in genetically modified crops new genes are introduced to the genome of the host in order to improve its resistance to certain pests/diseases or for improving the nutrient profile of the crop [42]. Another example is the commercial production of human proteins in recombinant microorganisms for therapeutic use [53,43,5,21]. This is sometimes based on the natural ability of certain bacteria to efficiently secrete properly folded human proteins (for example, insulin [21]). In this context, a fundamental problem is to maximize the translation rate of the heterologous gene (and thus the protein production rate) under the given constraints, e.g., the limited availability of intracellular components involved in translation. These constraints are needed also because very high initiation and elongation rates mean that the expression of the heterologous gene consumes too much resources of the translational machinery (e.g., ribosomes, tRNA molecules, etc), thus significantly deteriorating the fitness of the host. In addition, very high levels of protein abundance may eventually contribute to aggregation of proteins [45,31], leading to a decrease in the yield of heterologous protein production. All these aspects are encapsulated in the convex optimization problem that is addressed here for the RFM. We believe that this mathematical problem may thus be used to provide verifiable predictions on how to efficiently manipulate the various biological factors. There is a rich literature on using optimization theory, combined with evolutionary arguments, in biology (see, e.g., [54,46,2] and the references therein). This approach has been often criticized, but it has undoubtedly provided insight into the process of adaption under biological constraints, as well as helped to discriminate between alternative hypotheses for a suitable "fitness function" in various biological mechanisms. Furthermore, laboratory evolution experiments showed evolutionary adaptation of biological processes towards optimal operation levels. Examples include optimal metabolic fluxes in E. coli [26], and optimal protein expression levels from the lac operon [14]. We believe that the optimization problem posed here may lead to further progress in studying the evolution of the translation machinery.
The translation machinery is affected by mutations such as duplication/deletion of a tRNA gene or synonymous mutation affecting the codon bias usage. The concavity of the translation rate R may suggest that the selection of mutations that increase fitness indeed converges towards the optimal parameter values, as explained by the simple "hill climbing" argument described above (see Fig. 7).
Recent studies have shown that in various organisms the ribosomal density at the 5 ′ and 3 ′ ends of the ORF is higher than in the middle of the ORF (see, for example, [28,29,11,60]). In addition, the genomic ribosomal density is relatively constant in the middle of the ORF (usually more than 30 codons away from the two ORF ends). The elongation rate λ i is negatively correlated with the ribosomal density (or the probability that a site is occupied) at site i. Indeed, if λ i , that controls the elongation rate from site i, is small then there is a higher probability to see a ribosome in this site. Thus, these biological studies suggest that the elongation rates at the end of the chain are lower than in the middle of the chain, and that the rates near the middle are approximately equal. This agrees well with the optimal elongation rates derived based on our analysis in the case of equal weighting in the constraint (see Fig. 5).
Our results in the case of equal weights in the constraint also show that if the total biosynthetic budget b(n) is a sub-linear or linear function of n then R * decreases monotonically with n; however, when b(n) grows faster than a linear function in n then R * increases monotonically with n. The relation between expression levels and gene length has been studied experimentally. It has been shown that in some organisms, such as humans and S. cerevisiae [18,12], expression levels tend to monotonically decrease with gene length (shorter genes have higher expression levels). However, in other organisms, such as plants [50], an opposite relation was reported (longer genes have higher expression levels). Our analysis may suggest that one should take into account not only the difference in gene length, but also the difference in the available resources of the translational machinery.
An interesting question for further research is whether the translation rate in other models of translation, including various versions of TASEP [41,62,60,56,10,48], is also a concave function of its parameters.
Another possible research direction is the design and implementation of biological experiments based on the analytical results described above. Such experiments should combine: (1) methods for manipulating the translation machinery and/or the transcript of certain gene(s); and (2) online estimation of ribosomal density along the mRNA (e.g. using ribosome profiling [27]). The elongation rate of each codon can be estimated based on a method described in [13]. Manipulation of the translation machinery can include deletion of tRNA genes, and using the tet07 construct to down-regulate the initiation and elongation factors [20,6]. Techniques for local manipulation of a transcript include generating libraries of a certain heterologous non-functional gene (e.g., a GFP protein). In each of the variants a few mutations (relatively to the wild-type) are introduced either in the 5'UTR (corresponding to λ 0 ) or the ORF (corresponding to λ 1 , λ 2 , . . . ), and the protein levels and ribosomal densities are measured [33,65]. The fact that the heterologous gene is non-functional to the host assures that the observed changes in translation efficiency are due to the introduced modifications.
As a specific example, one can measure the effect of modifying the elongation rates of different codons (corresponding to λ 1 , λ 2 , ..) by introducing synonymous mutations in different parts of the ORF. We expect that a graph depicting the translation rate as a function of elongation rates will be concave (as in Fig. 4).
Finally, the effect of single mutations in different parts of the transcript on translation rate is a fundamental question related to various biomedical disciplines. Specifically, it is known that codon substitutions in different parts of the coding sequence affect elongation and initiation rates (i.e. the λ i s) via various mechanisms (e.g. mRNA folding and adaptation to the tRNA pool [16,22,33,62]). Our result is based on linking R to the Perron root of a tridiagonal matrix that depends on the λ i s. This can serve as a starting point for sensitivity analysis of R, i.e. analyzing the effect of small changes in the λ i s on R. This topic is currently under study.
Step 3: Recall that the spectral radius of a square matrix is the maximum over the absolute values of its eigenvalue. The spectral radius of a non-negative matrix is an eigenvalue of the matrix called the Perron root [19]. The next result shows that R −1/2 is the largest eigenvalue of the non-negative matrix A.
Proof of Proposition 1. Let v = v 1 . . . v n+2 ′ ∈ R n+2 + denote a Perron eigenvector of the symmetric matrix A, i.e., an eigenvector corresponding to the Perron root R −1/2 . It follows from known results (see, e.g., [37]) that and combining this with (26) yields Since all the components of the Perron eigenvector are strictly positive [25], this implies that ∂ ∂λi R > 0 for all i = 0, ..., n. Proof of Theorem 2. The proof is based on formulating the Lagrangian function associated with Problem 1 and determining the optimal parameter values by equating its derivatives to zero (see, e.g., [8]). The Lagrangian is L(λ 0 , . . . , λ n , θ) := R(λ 0 , . . . , λ n ) + θ b − n i=0 w i λ i , where θ ∈ R is the Lagrange multiplier. Differentiating this with respect to λ i and equating to zero yields where | * means that the equation holds once the optimal values λ * i are substituted. This implies in particular that and combining this with (30) yields where v * is the unique (up to scaling) Perron eigenvector of the non-negative and irreducible matrix A * := A(λ * 0 , . . . , λ * n ) [25,Ch. 8]. The equation A * v * = (R * ) −1/2 v * yields , and substituting this in (31) and simplifying yields (13).