Deciphering mRNA Sequence Determinants of Protein Production Rate
Juraj Szavits-Nossan Luca Ciandrini M. Carmen Romano

	



jszavits@staffmail.ed.ac.uk
SUPA, School of Physics and Astronomy, University of Edinburgh, Peter Guthrie Tait Road, Edinburgh EH9 3FD, United Kingdom


DIMNP UMR 5235, Université de Montpellier and CNRS, F-34095, Montpellier, France
Laboratoire Charles Coulomb UMR5221, Université de Montpellier and CNRS, F-34095, Montpellier, France


SUPA, Institute for Complex Systems and Mathematical Biology, Department of Physics,  Aberdeen AB24 3UE, United Kingdom
Institute of Medical Sciences, University of Aberdeen, Foresterhill, Aberdeen AB24 3FX, United Kingdom


One of the greatest challenges in biophysical models of translation is to identify coding sequences features that affect the rate of translation and therefore the overall protein production in the cell. We propose an analytic method to solve a translation model based on the inhomogeneous totally asymmetric simple exclusion process, which allows us to unveil simple design principles of nucleotide sequences determining protein production rates. Our solution shows an excellent agreement when compared to numerical genome-wide simulations of S. cerevisiae transcript sequences and predicts that the first 10 codons, which is the ribosome footprint length on the mRNA, together with the value of the initiation rate, are the main determinants of protein production rate under physiological conditions. Finally, we interpret the obtained analytic results based on the evolutionary role of codons' choice for regulating translation rates and ribosome densities.


87.16.aj, 87.10.Mn, 05.60.-k


















Deciphering mRNA Sequence Determinants of Protein Production Rate
    M. Carmen Romano
    
=================================================================







Translation is one of the major steps in protein biosynthesis. During this process, the nucleotide sequence of a messenger RNA (mRNA) is translated into a functional protein. Each nucleotide triplet, called codon, codes for a specific amino acid, the proteins' building block. There is experimental evidence that the rate at which a certain mRNA is translated depends on its specific codon sequence <cit.>,  especially in the case of eukaryotes. Identifying sequence features that determine protein production rate, also commonly referred to as translation rate or efficiency, is a fundamental open question in molecular biology <cit.>.
	

	
Translation is performed by molecular motors called ribosomes, which move unidirectionally along the mRNA. The amino acids are delivered to the ribosome by molecules called transfer RNAs (tRNAs), which are specific to the codon and the amino acid they deliver (Fig. <ref>). The dwelling time of a ribosome on a specific codon depends primarily on the abundance of the corresponding tRNA <cit.>. Understanding how codon sequences determine protein production rates can potentially unlock many synthetic applications <cit.>.
	
The standard biophysical model of translation is known as the totally asymmetric simple exclusion process (TASEP), which captures the concurrent motion of ribosomes on the mRNA <cit.>. In this model ribosomes progress along the mRNA codon by codon, provided that the codon a ribosome moves onto is not occupied by another ribosome. The protein production rate can then be identified as the ribosomal current of this driven lattice gas. Due to the net current, TASEP is not in equilibrium and its steady state is only known for a few special cases <cit.>. Unfortunately, most biologically relevant variants of TASEP can be studied only numerically <cit.>, and efficient methods of exploring a large number of parameters are lacking. The steady state of the TASEP with non-uniform hopping rates is a long outstanding problem in nonequilibrium statistical physics <cit.>. An analytic prediction of protein production rates is also needed in order to interpret recently developed ribosome profiling experiments that are capable of monitoring ribosome positions along the mRNA <cit.>. Such analytical prediction is missing in previous models of ribosome dynamics implicitly or explicitly based on the TASEP <cit.>, and it is key to decipher sequence determinants of protein production rates.
	
In this Letter, we develop a versatile analytic method to solve TASEP-based models of translation. Our analytic approach, integrated with simulations and experimental data, allows us to efficiently identify the main features of mRNA codon sequence that determine the rate of protein production. The genome-wide comparison of our analytic predictions with numerical simulations shows an excellent agreement for the model organism S. cerevisiae (baker's yeast).  
	




Stochastic model of mRNA translation. We focus on the model for translation introduced in <cit.>. The mRNA is represented by a one-dimensional lattice consisting of L discrete sites (codons), where site 1 designates the start codon. Ribosomes are represented by particles that occupy ℓ = 10 lattice sites, which is the ribosome footprint length measured in ribosome profiling experiments  <cit.>. We identify the position 1≤ i≤ L of a ribosome with the position of its A-site, which is located d = 5 lattice sites from the trailing end of the ribosome (Fig. <ref>a). A ribosome that waits for a tRNA at position i is labelled by 1_i and a ribosome that has already received the correct tRNA and is ready to move is labelled by 2_i. The set of labels of all translating ribosomes on the lattice is called a configuration C of the system. For example, C=1_1 2_13 denotes a lattice with 2 ribosomes, one at site 1 waiting for a tRNA and another one at site 13 that has already received the correct tRNA.

Ribosomes initiate translation at rate α by binding to the mRNA so that their A-site is at the start codon, provided the sites 1,…,ℓ-d+1 are empty (Fig. <ref>a). A ribosome at site i makes the transition 1_i → 2_i at rate k_i dependent on tRNA abundances (Fig. <ref>b) and moves one site forward at rate γ (Fig. <ref>c). Due to steric interactions between particles, the particle at site i can move only if there is no ribosome at site i+ℓ. Termination occurs when the particle at site L receives the last amino acid, releases the final protein and detaches from the mRNA, which we integrate into a single step occurring at rate β (Fig. <ref>d). 

Our main goal is to compute the rate of protein production as a function of the parameters of the model, which are α, β, γ and  k_i for each of the L codons used. We assume that translation takes place under steady-state conditions, so that the ribosomal current is constant along the mRNA and is equal to the rate of protein production, which in turn is equal to rate at which ribosomes load onto the mRNA and initiate translation. To this end we define the codon occupation number τ_i to be equal to 1 if the i-th codon is occupied by an A-site, and 0 if otherwise. By definition, the exact steady-state ribosomal current J then reads


    J=α∑_C[∏_i=1^ℓ(1-τ_i(C))]P(C),


where τ_i(C) denotes the i-th codon occupation number for configuration C, P(C) denotes the steady-state probability that the lattice is in configuration C, and the summation goes over all configurations C. Other quantities of interest that we compute are the local and total particle densities ρ_i=⟨τ_i⟩ and ρ=(1/L)∑_i=1^Lρ_i, respectively.





Series expansion method for computing P(C). In order to find P(C) one has to solve the steady-state master equation M𝐏=0, where 𝐏 is a column vector whose 𝒩 elements are the steady-state probabilities P(C) of being in configuration C, and 𝒩 denotes the total number of configurations. The transition rate matrix M_C,C' is given by W_C'→ C for C≠ C' and -e(C) for C=C', where W_C'→ C is the transition rate from C' to C and e(C)=∑_C”≠ CW_C→ C” is the total exit rate from C. The exact solution of the master equation can be formally written as


    P(C)=detM^(p,p)/∑_q=1^𝒩detM^(q,q),


where p is the position of configuration C in the column vector 𝐏 and detM^(p,p) is a determinant of the matrix obtained by removing p-th row and p-th column from M (see Supplemental Material for details). Unfortunately, calculating this determinant is feasible only for unrealistically small system sizes.

To circumvent this problem, we exploit the fact that detM^(p,p) is a multivariate polynomial in the variables α, k_1,…,k_L, β and γ <cit.>. In the biological literature it is often assumed that the initiation rate α is a major limiting step of the translation process, mainly determined by the presence of secondary structures <cit.>. For this reason we use α as an expansion parameter and we assume that α≪γ,β,k_1,…,k_L. By collecting terms with the same power of α, we can rewrite P(C) as an univariate polynomial f(C) in the variable α with unknown coefficients f_n(C)


    P(C)=f(C)/∑_Cf(C),   f(C)=∑_n=0^K(C)f_n(C)α^n,


where the coefficients f_n(C) depend on the transition rates k_1,…,k_L,γ and β (in order to ease the notation, we leave out the explicit dependence on those parameters). Since we expect that f(C) can be well approximated by the first few terms, the value of K(C) in Eq. (<ref> is irrelevant in our study.

In order to find the unknown coefficients f_n(C), we insert Eq. (<ref>) into the master equation M𝐏=0, collect all the terms with the same n-th power of α and equate their sum to zero. For a given power n, the resulting equation is similar to the original master equation in which P(C) is replaced by f_n(C) unless the coefficient multiplying P(C) is α, in which case P(C) is replaced by f_n-1(C). Starting with n=0, we note the equations for the coefficients f_0(C) have the same form as the original master equation, but with α=0, leading to the trivial solution f_0(C)=0 if C≠∅. Since the terms f(C) are not normalized, we have the freedom to choose any value for f_0(∅), which we set to 1.

For n≥ 1, the equation for f_n(C) reads


    e_0(C)f_n(C)+f_n-1(C)∑_C'I_C,C'=∑_C'I_C',Cf_n-1(C')
         +∑_C'(1-I_C',C)W_C'→ Cf_n(C') ,


where I_C,C'=1 if W_C→ C'=α and is 0 otherwise, and  e_0(C)=∑_C'(1-I_C,C')W_C→ C' is the total exit rate from C excluding the rate α. 

A key observation in our analysis is that f_n(C)=0 whenever the number of particles in C is larger than n. This follows from the result f_0(C)= 1   (0) if C=∅ (C≠∅) in conjunction with the hierarchical structure of Eq. (<ref>), which connects configurations differing in the number of particles by no more than one (a proof for n=1 is presented in the Supplemental Material). This allows us to write Eq. (<ref>) taking into account configurations with only one particle and discarding all the others, which yields


	
	
    f_1(1_1)=1/k_1f_0(∅),   f_1(1_L)=β/k_Lf_1(2_L)
        f_1(1_i)=γ/k_i f_1(2_i-1),   2 ≤ i≤ L,
        f_1(2_i)=k_i/γ f_1(1_i),   1≤ i≤ L-1.



The solution to Eqs. (<ref>) is given by


    f_1(1_i)=1/k_i,   f_1(2_i)=1/γ    i=1,…,L-1
    1/β    i=L


The equations for f_2(C) involving configurations with one and two particles are presented in the Supplemental Material. 

Once we determine the coefficients f_n(C) up to a desired order n, we can compute the steady-state average of any observable 𝒪(C) by inserting P(C) from Eq. (<ref>) and expanding ⟨𝒪⟩ around α=0,


    ⟨𝒪(C)⟩=∑_n=0^K∑_C𝒪(C)f_n(C)α^n/∑_n=0^K∑_Cf_n(C) α^n=∑_n=0^∞c_nα^n.


For example, the first three coefficients c_0, c_1 and c_2 are given by


    c_0=a_0/b_0,  c_1=a_1-c_0 b_1/b_0,  c_2=a_2-c_1 b_1-c_0 b_2/b_0,


where a_n and b_n are defined as


    a_n=∑_C𝒪(C)f_n(C),    b_n=∑_Cf_n(C).


Notice that the expansion in Eq. (<ref>) is slightly different for the current J due to an extra α in Eq. (<ref>) and is given by J=∑_n=0^∞c_n α^n+1. Using the expressions for f_n(C) and f_n(C) computed earlier yields 


    J=α[1-∑_i=1^ℓ(1/k_i+1/γ_i)α+O(α^2)],
       ρ=1/L∑_i=1^L(1/k_i+1/γ_i)α+O(α^2),
       ρ_i= (1/k_i+1/γ_i)α+O(α^2),


where γ_i=γ+(β-γ)δ_i,L. These equations constitute our main result. Equation (<ref>) shows that if the initiation rate α is small compared to k_1,…,k_L and γ, than the protein production rate J depends predominately on the initiation rate, along with the translocation and elongation rates of the first 10 codons, corresponding to the ribosome footprint ℓ. The importance of the first 10 codons is a direct consequence of the excluded volume interactions: any ribosome already present in that region will prevent a new ribosome from binding the mRNA. It is also important to emphasize that (<ref>)–(<ref>) are exact series expansions around α=0; the approximation is made only when the series is truncated.

Independent Particle Approximation (IPA). Interestingly, the excluded volume interactions between particles have no effect on the first two terms in the series expansion of J. This motivates us to ask how the expansion in Eq. (<ref>) would look like if we assumed that all particles are independent, i.e. not experiencing any exclusion interaction. In our model, the IPA amounts to replacing P(C) with


    P^IPA(C)=1/Z_L∏_j=1^N(C)w_1(θ_X(j)),


where N(C) is the number of particles in a configuration C, θ is one of the two particle states 1 and 2, X(i) is the position of the i-th particle on the lattice and Z_L is the normalization constant. The weights w_1(1_i) and w_1(2_i) for i=1,…,L are obtained by solving the master equation for a single particle and are given by α/k_i and α/γ_i, respectively. The corresponding expressions for the current J and local density ρ_i read


    J^IPA=α/∏_i=1^ℓ(1+p_i),  ρ^IPA_i=p_i/1+p_i,


where p_i=α(1/k_i+1/γ_i). We will use these results later in order to determine the importance of ribosome collisions in real genetic sequences. We also note that the IPA in (<ref>) provides a good approximation to f_n(C) for n=N(C) when the particles in C are far apart from each other.





Application to mRNA translation in yeast. We now apply our results to the transcriptome of S. cerevisiae using realistic model parameters. The values of α in the range 0.005-4.2 s^-1 with the median value of 0.09 s^-1 have been previously estimated in Ref. <cit.> using genome-wide experimental values of the ribosomal density. We assume that the rates k_i are mainly proportional to the gene copy number of tRNAs delivering the corresponding amino acid <cit.>; the rates are normalized so that the average codon translation rate is equal to the experimental value of 10 codons/s <cit.>; the estimates of all elongation rates along with the distribution of α is presented in the Supplemental Material. The translocation rate γ is fixed to γ=35 codons/s <cit.>, and termination is assumed to be fast and comparable to translocation <cit.>, β≈γ, so that γ_i=γ ∀ i in Eqs. (<ref>)-(<ref>). 

In total, we analyzed 5836 gene sequences; for each gene we calculated J, ρ and ρ_i for 1≤ i≤ L up to and including the second order of the perturbative expansion at the corresponding physiological value of α. The results were then compared to the exact values obtained numerically with stochastic simulations using the Gillespie algorithm <cit.> by calculating the percent error ϵ.



For the protein production rate J, with the zeroth order of the perturbative expansion (predicting that J=α) we obtain an error of ϵ < 5% for only 11% of the genes. Remarkably, that percentage jumps to 80.7% when the first-order coefficients in Eq. (<ref>) are taken into account (Fig. <ref>). Including the second-order coefficients (computed numerically from Eq. (<ref>)) does not significantly improve results, due to a large value of α⪆ 0.15 s^-1 in about 20% of the genes. Since the coefficients in Eq. (<ref>) typically alternate in sign, truncating the series will ultimately lead to a wrong result when the value of α is large enough. On the other hand, the IPA does not suffer from this problem and leads to an error ϵ<5% in 94% of genes, whereby only 1% of genes have ϵ>20%. The success of the IPA also suggests that ribosome collisions and traffic jams have a minor effect on the rate of translation, which is in accordance with recent experimental evidence <cit.>. This is also apparent from the density profile ρ_i, which is very well approximated by the linear approximation in Eq. (<ref>), even for larger values of α (Fig. <ref> b and c).





Identifying determinants of mRNA translation. Our analytic prediction allows us to decompose the contributions from initiation and elongation to the rate of translation J, thereby addressing a long-standing question about main determinants of protein production rate. Remarkably, the expressions for the current obtained with both the first-order and the independent particle approximation involve only the first ℓ=10 codons. This result therefore strongly indicates that, together with the initiation rate, the first 10 codons of the mRNA are the key determinants of the protein production rate by preventing a new ribosome from binding which effectively decreases the initiation rate.



If we assume that the cell maximizes the rate of protein production, given the above result, we would expect to find a signature in the genome for selecting efficient fast codons at the beginning of each gene. We test this hypothesis by computing  the ribosomal current for the fastest (J^F) and the slowest (J^S) set of first ℓ synonymous codons for each gene, using the IPA. A score η_J=(J-J^S)/(J^F-J^S) is then assigned to each gene, which is 1 (0) when the sequence corresponds to the fastest (slowest) codon sequence. On the other hand, one might assume that the cell not only tries to maximize protein production rates, but at the same time it tries to minimize the ribosome density ρ on mRNAs. This assumption is motivated by the fact that ribosomes are limiting <cit.> and highly costly in terms of cellular energy resources <cit.>, and therefore ribosome queues are to be avoided. Hence, we also compute η_ρ=(ρ-ρ^S)/(ρ^F-ρ^S) for each gene, where ρ^F and ρ^S denote the ribosome density for the fastest and slowest set of synonymous codons, respectively.

Figure <ref> shows the histogram of η_J (yellow) and η_ρ (blue) computed for 5836 genes of S.cerevisiae. Both histograms show an average of 0.7 suggesting the selection of fast codons near the start codon to maximize J, as well as an overall selection of fast codons along the mRNA to minimize ρ. However, the width of the distribution of η_ρ is substantially smaller than the one of η_J. This might indicate that the optimization of protein production rate is strongly dependent on the particular gene, since different proteins are needed at different concentrations. In contrast, the minimization of the number of ribosomes on mRNAs could be a more general constraint. 





Conclusions. We have presented an analytic method that allows us to quantify sequence determinants of protein production rates, using a model for translation that is based on an inhomogeneous exclusion process. Our results demonstrate that the rate of protein production is largely determined by the initiation rate and the elongation rates of the first 10 codons (the ribosome footprint length on the mRNA), which control how fast ribosomes load onto the mRNA, whereby ribosome collisions and queues have a minor effect under physiological conditions.






Acknowledgments. MCR and LC contributed equally to this work. JSN was supported by the Leverhulme Trust Early Career Fellowship. MCR was supported by the Biotechnology and Biological Sciences Research Council (BBSRC) BB/N017161/1 and the Scottish Universities Life Sciences Alliance. LC thanks the CNRS for having granted him a “demi-dé́légation” (2017–18).






99
	
Gingold11 H. Gingold, and Y. Pilpel,  Mol. Syst. Biol. 7, 481 (2011).
Kemp13 A.J. Kemp, R. Betney, L. Ciandrini, A. C. Schwenger, M. C. Romano, and  I. Stansfield, Mol. Microbiol. 87(2), 284 (2013).
Chu14 D. Chu, E. Kazana, N. Bellanger, T. Singh, M. F. Tuite, and T. von der Haar, EMBO J. 33(1), 21 (2014).
Gorgoni16	B. Gorgoni, L. Ciandrini, M. R. McFarland, M. C. Romano, and I. Stansfield, Nucleic Acids Res. 44(19), 9231 (2016).
Varenne84 S. Varenne, J. Buc, R. Lloubes, and C. Lazdunski, J. Mol Biol. 180(3), 549 (1984).
Sorensen89 M. A. Sorensen, C. G. Kurland, and S. Pedersen, J. Mol. Biol. 207, 365–377 (1989).
Gorochowski16 T. E. Gorochowski, I. Avcilar-Kucukgoze, R. A. Bovenberg, J.A. Roubos, and Z. Ignatova,  ACS Synth. Biol. 5(7), 710 (2016).
Boel16 G. Boël et al, Nature 529(7586), 358 (2016).
MacDonald68 C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin, Biopolymers 6, 1–25 (1968); C. T. MacDonald and J. H. Gibbs, Biopolymers 7, 707–25 (1969).
Tuller16 H. Zur, and T. Tuller, Nucleic Acids Res. 44(19), 9031–9049 (2016).
DEHP B. Derrida, M. R. Evans, V. Hakim, and V. Pasquier, J. Phys. A: Math. Gen. 26, 1493 (1993); G. M. Schütz, and E. Domany, J. Stat. Phys. 72, 277–96 (1993).
BlytheEvans07 R. A. Blythe, and M. R. Evans, J. Phys. A: Math. Theor. 40, R333–R441 (2007).
ChouLakatos04 T. Chou, and G. Lakatos, Phys. Rev. Lett. 93, 198101 (2004).
Shaw04 L. B. Shaw, J. P. Sethna, and K. H. Lee, Phys. Rev. E 70, 021901 (2004).
Zia11 R. K. P. Zia, J. J. Dong, and B. Schmittmann, J. Stat. Phys. 144, 405 (2011).
Schmidt15 J. Schmidt, V. Popkov and A. Schadschneider, EPL 110, 20008 (2015)
Ingolia09  N. T. Ingolia, S. Ghaemmaghami†, J. R. S. Newman, and J. S. Weissman, Science 324, 218–223 (2009).
Gilchrist06 M. A. Gilchrist, A. Wagner, J. Theor. Biol., 239(4), 417-434 (2006).
Mitarai08 N. Mitarai, K. Sneppen, and S. Pedersen, J. Mol. Biol., 382(1), 236-245 (2008).
Reuveni11 S. Reuveni, I. Meilijson, M. Kupiec, E. Ruppin, and T. Tuller, PLoS Comput. Biol. 7, e1002127 (2011).
Ciandrini13 L. Ciandrini, I. Stansfield, and M. C. Romano, PLoS Comput. Biol. 9(1), e1002866 (2013).
Shah13 P. Shah, Y. Ding, M. Niemczyk, G. Kudla, and J.B. Plotkin, Cell 153(7), 1589 (2013).
Ciandrini10 L. Ciandrini, I. Stansfield, and M. C. Romano. Phys. Rev. E 81, 051904 (2010).
Klumpp08 S. Klumpp, Y. Chai, and R. Lipowsky, Phys. Rev. E 78, 041909 (2008).
SzavitsNossan13  J. Szavits-Nossan, J. Phys. A: Math. Theor. 46, 315001 (2013).
Kudla09 G. Kudla, A. W. Murray, D. Tollervey, and J. B. Plotkin, Science 324, 255–258 (2009).
Salis09 H. M. Salis, E. A. Mirsky, and C. A. Voigt, Nature Biotechnol. 27, 946–950 (2009).


Percudani97 R. Percudani, A. Pavesi, and S. Ottonello, J. Mol. Biol. 268(2), 322–330 (1997).


Savelsbergh03 A. Savelsbergh, V. I. Katunin, D. Mohr, F. Peske, M. V. Rodnina, and W. Wintermeyer, Mol. Cell 11, 1517–1523 (2003).


Arava03 Y. Arava, Y. Wang, J. D. Storey, C. L. Liu, P. O. Brown, and D. Herschlag, PNAS 100(7), 3889–3894 (2003).


Gillespie77 D. T. Gillespie, J. Phys. Chem. 81(25), 2340 (1977).

Laursen10 B. S. Laursen, H. P. Sørensen, K. K. Mortensen, and H. U. Sperling-Petersen, Microbiol. Mol. Biol. Rev. 69, 101 (2005).
Greulich12 P. Greulich, L. Ciandrini, R. J. Allen, and M. C. Romano. Phys. Rev. E 85, 011142 (2012).
Keiler15 K. C. Keiler, Nature Rev. Microbiol. 13, 285 (2015).
	












	Supplemental Material to:
 Deciphering mRNA Sequence Determinants of 
 Protein Production Rate




[1]













§ EXACT SOLUTION OF THE STEADY-STATE MASTER EQUATION


The exact solution, Eq. (3) in the main text, pertains to any ergodic Markov jump process with a finite number of states. In order to derive the exact solution, let us write again the steady-state master equation in a matrix form


    M𝐏=0,


where 𝐏 is a column vector whose elements are steady-state probabilities P_i, the 𝒩 states are indexed by i=1…,𝒩 and the 𝒩×𝒩 matrix M is given by


    M_ij=
    		W_jk    i≠ j
    
    		-∑_k≠ iW_ik    i=j,


where W_ij is a transition rate from a state i to j. We note from Eq. (<ref>) that the sum of all elements in each column of M is zero


    ∑_i=1^𝒩M_ij=∑_i=1  i≠ j^𝒩W_ji-∑_k=1  k≠ j^𝒩W_jk=0.


Consequently, the sum all row vectors of M is zero, which means that the vectors are linearly dependent and thus the determinant of M is equal to zero,  


    detM=0.
 

Combining Eq. (<ref>) with the Laplace expansion of a determinant yields


    0=detM=∑_i=1^𝒩M_ijC_ij=∑_j=1^𝒩M_ijC_ij.


Here C_ij is a cofactor of M defined as C_ij=(-1)^i+jdetM^(i,j), where the matrix M^(i,j) is obtained from M by removing i-th row and j-th column . Inserting Eq. (<ref>) into (<ref>) gives


    detM   =∑_i=1  i≠ j^𝒩M_ijC_ij+M_jjC_jj=∑_i=1  i≠ j^𝒩M_ijC_ij-∑_k=1  k≠ j^𝒩M_kjC_jj
       =∑_i=1  i≠ j^𝒩M_ij(C_ij-C_jj)=0,


from which we conclude that 


    C_ij=C_jj,


for any i and j, i.e. the cofactor C_ij does not depend on the state i. Inserting Eq. (<ref>) back into Eq. (<ref>) gives


    ∑_j=1^𝒩M_ijC_jj=0,


which is precisely the starting steady-state master equation, Eq. (<ref>). We thus conclude that


    P_i=C_ii/∑_jC_jj=detM^(i,i)/∑_j=1^𝒩detM^(j,j).




§ ARGUMENT FOR F^(N)(C)=0 WHENEVER THE NUMBER OF PARTICLES IN A CONFIGURATION C IS LARGER THAN N


As noted in the main text, a key observation in our analysis is that f^(n)(C)=0 whenever the number of particles in C is larger than n. We present here the argument for the case of n=1; a similar argument applies to higher orders. For n=1, using Eq. (5) in the main text we obtain


    f^(1)(C)=1/e(C)[∑_C'I_C'→ Cδ_C',∅+∑_C'(1-I_C'→ C)W_C'→ Cf^(1)(C') - δ_C,∅].


Consider C to be one full lattice configuration with all N=⌊ (L-1+ℓ)/ℓ⌋ particles in state 1. The first term on the right hand side of Eq. (<ref>) is zero because there is no direct transition from the empty configuration ∅ to C. The second term on the right hand side is also zero, because C can only be accessed through the transition rate α. Finally, the third term on the right hand side is clearly zero. Hence, f^(1)(C)=0. Now, it is easily seen that all full lattice configurations C' that can be accessed only from C will fulfil f^(1)(C')=0. Iterating further this procedure to the configurations that can be accessed from those, it becomes clear that all full lattice configurations C have f^(1)(C)=0. This argument can be further extended to all configurations with N-1, N-2, N-3, ..., 2 particles, resulting in f^(1)(C)=0 for all of them. 

The situation changes when we consider configurations with only 1 particle on the lattice, since then there is a direct transition from the empty configuration ∅ to one of those configurations, leading to f^(1)(C) ≠ 0. The rest of the coefficients f^(1)(C) for 1-particle configurations will depend on each other, and hence, they do not vanish. Therefore, we can conclude that only 1-particle configurations have non-vanishing f^(1)(C). For n=1, this key observation allows us to write Eq. (5) taking into account configurations with only one particle and discarding all the others, which yields Eqs. (6a)-(6c) in the main text.



§ EQUATIONS FOR THE SECOND-ORDER COEFFICIENTS F^(2)(C)


As stated in the main text, all second-order coefficients f^(2)(C) whereby a configuration C has more than two particles are equal to zero. The equations for the non-zero coefficients f^(2)(C) are presented below, where we use the the notation γ_j=γ for 1≤ j≤ L-1 and γ_L=β.

First, we look at configurations with one particle at site i=1 and the other at site j=i+l,…,L:


	
	
    f^(2)(1_1 1_j)=1/k_1+k_j[f^(1)(1_j)+(1-δ_1+l,j)γ_j-1f^(2)(1_1 2_j-1)]
        f^(2)(1_1 2_j)=1/k_1+γ_j[f^(1)(2_j)+k_j f^(2)(1_1 1_j)]
        f^(2)(2_1 1_j)=1/(1-δ_1+l,j)γ_1+k_j[k_1 f^(2)(1_1 1_j)+(1-δ_1+l,j)γ_j-1f^(2)(2_1 2_j-1)]
        f^(2)(2_1 2_j)=1/(1-δ_1+l,j)γ_1+γ_j[k_1 f^(2)(1_1 2_j)+k_j f^(2)(2_1 1_j)].



where we used the Kronecker delta to account for the excluded volume interaction. Next, we look at configurations with particles at sites i=2,…,L-l and j=i+l,…,L.


	
	
    f^(2)(1_i 1_j)=1/k_i+k_j[γ_i-1f^(2)(2_i-1 1_j)+(1-δ_i+l,j)γ_j-1f^(2)(1_i 2_j-1)]
        f^(2)(1_i 2_j)=1/k_i+γ_j[γ_i-1f^(2)(2_i-12_j)+k_j f^(2)(1_i 1_j)]
        f^(2)(2_i 1_j)=1/(1-δ_i+l,j)γ_i+k_j[k_i f^(2)(1_i 1_j)+(1-δ_i+l,j)γ_j-1f^(2)(2_i 2_j-1)]
        f^(2)(2_i 2_j)=1/(1-δ_i+l,j)γ_i+γ_j[k_i f^(2)(1_i 2_j)+k_j f^(2)(2_i 1_j)].



Finally, the equations for configurations with only one particle are given by


	
	
    f^(2)(1_1)=1/k_1[γ_L f^(2)(1_1 2_L)+f^(1)(∅)],
        f^(2)(2_1)=k_1/γ_1f^(2)(1_1)+γ_L/γ_1f^(2)(2_1 2_L),
        f^(2)(1_i)=1/k_i[γ_i-1f^(2)(2_i-1)+θ(L-l+1-i)γ_L f^(2)(1_i 2_L)-f^(1)(1_i)],   i=2,…,L
        f^(2)(2_i)=1/γ_i[k_i f^(2)(1_i)+θ(L-l+1-i)γ_L f^(2)(2_i 2_L)-f^(1)(2_i)],   i=2,…,L



where θ(n)=0 for n<0 and 1 for n≥ 0. 

The equations (<ref>)-(<ref>) can be easily solved numerically by iteration. We first solve Eqs. (<ref>) for i=1,j=l+1. We then iterate the recursion relation in Eqs. (<ref>) for fixed i=1 and j=2+l,…,L. We then solve equations (<ref>) for fixed i=2 and j=l+3,…,L and repeat this procedure until i=L-l and j=L. Finally, we calculate the coefficients f^(2)(1_i) and f^(2)(2_i) for the single-particle configurations using f^(2)(1_i 2_L) and f^(2)(2_i 2_L) that we solved in the previous steps. We note that the only unknown that we cannot determine is f^(1)(∅), but it turns out that all terms containing f^(1)(∅) will cancel out later when we calculate the coefficients in the series expansions (11)-(13) in the main text.



§ ESTIMATES OF THE CODON ELONGATION RATES


The total translation rate ω_i of the codon i is:


    1/k_i+1/γ=1/ω_i,


where the translocation rate γ is set to 35s^-1 as mentioned in the main text. The codon elongation rate k_i, which represents the average arrival and recognition time of the cognate tRNA, is determined by following the procedure introduced in <cit.>, which we report here. 

For each of the 41 tRNAs types j we consider their gene copy number (GCN), which allows us to provide a first estimate of the rate k_j:


    k_j=rGCN_j/∑_j=1^41 GCN_j,


where GCN_j is the gene copy number of the tRNA of type j with j=1,…,41, and r is a proportionality constant. These rates were then adjusted to take into account experimental evidence suggesting that the translation rates of codons using the G-U wobble are reduced by 39% compared to their G-C counterparts; analogously, codons using the wobble I-C and codons using the wobble I-A are reduced by 36% relative to their I-U counterparts <cit.>. To calculate the proportionality constant r, we used the experimental value of 10 codons/s for the average codon translation rate ⟨ω_i ⟩ defined as


    ⟨ω_i ⟩=∑_i=1^61( k_i γ/k_i+γ)n_i/n


where n_i/n is the relative abundance of all different codon types in the cell, n_i is total number of codons in the cell of exactly type i, and n=∑_i=1^61n_i (the three STOP codons are excluded) <cit.>. Table <ref> summarises the resulting elongation rates k_i for 61 codons. 



The estimates provided are valid in physiological conditions, when amino acids are not limiting. In case of amino acid starvation the rates of the codons affected should be modified (the elongation rates depend on the abundance of tRNAs charged with the correct amino acid).



§ DISTRIBUTION OF THE TRANSLATION INITIATION RATES


The translation initiation rates α were estimated for each gene in Ref. <cit.> by matching the average ribosome density predicted by the model with the experimental value of the ribosome density obtained in Ref. <cit.>. The resulting values for 5836 genes that we analyzed in the main text are in the range 0.005-4.2 s^-1 with the median value of 0.09 s^-1. The histogram of the rates is presented in Figure <ref>, showing that the initiation rate α for the majority of genes is indeed much smaller than the smallest elongation rate, which has the value of 1.02 codons/s for the CUU codon (see Table <ref>). The list of all translation initiation rates can be found in the Supplemental Material of Ref. <cit.>.




10

Ciandrini13 L. Ciandrini, I. Stansfield, and M. C. Romano, PLoS Comput. Biol. 9(1), e1002866 (2013).
	
Gilchrist06
	M. A. Gilchrist and A. Wagner, J. Theor. Biol. 239 417–434 (2006).
	
Arava03 Y. Arava, Y. Wang, J. D. Storey, C. L. Liu, P. O. Brown, and D. Herschlag, PNAS 100(7), 3889–3894 (2003).
	



