Stochastic model of ribosome drop-off

Review of TASEP.

The original TASEP model [36] consists of a one-dimensional lattice of N sites spanning a length L. Particles hop stochastically along this lattice, from site to site. The particles have hard sphere properties, and therefore can neither occupy the same site, nor can they overtake each other. The dynamics are defined through the microscopic behaviour of particles, which can

1. enter at a rate α, if site 1 is unoccupied (left boundary),
2. hop from site i to site i + 1 at a rate k, if site i + 1 is empty (bulk),
3. exit the lattice at a rate β from site N (right boundary).

The language above, used commonly, is adapted to a transport problem in which a current of particles flows from the left (upstream) boundary to the right (downstream) boundary. In the context of translation, the lattice sites correspond to codons of the mRNA strand, along which ribosomes advance stochastically as they translate the sequence, from the left (5’) end to the right (3’) end. The microscopic rates correspond to initiation (α), elongation (k) and termination (β).

The state of the system is parameterised in terms of its occupation numbers n_i(t) for each site i of the lattice, so that n_i(t) is either 1 (occupied) or 0 (empty) at time t. We can then write master equations for the time-averaged occupation numbers ρ_i = 〈n_i(t)〉_t, which in the mean-field approximation (neglecting correlations between the sites) is given by (1) (2) (3) (Note that the Ribosome Flow Model (RFM) [40] analyses the solutions of the mean-field approximation of TASEP-like models of translation.) We shall assume a steady state for all ρ_i. We obtain the current of particles at a given position (number of particles going from site i to site i + 1 per unit time) as (4) As the number of particles is conserved along the lattice, this current does not depend on the site index i. Notice that every time a particle (ribosome) advances by one site (codon), another element is added to the nascent chain of amino acids. Thus the steady-state current is proportional to the amino acid incorporation rate and, in particular, the current at the 3’ end (right hand side) is proportional to the number of fully formed proteins produced per unit time, i.e. to the translation rate.

Four different phases [42–44] can be distinguished, depending on the relative values of α, β and k, and characterised by the density of particles along the lattice:

• low density (LD): if α < β and α < k/2, the entry rate limits the dynamics. The density in the bulk is ρ = α/k and J = α(1 − α/k).
• high density (HD): if β < α and β < k/2, the exit rate is the limiting step of the process. The density in the bulk is given by ρ = 1 − β/k and J = β(1 − β/k).
• maximal current (MC): if α, β > k/2 the system is limited by the bulk hopping rate k, and it achieves the maximal possible current J = k/4 at a density ρ = 1/2.
• shock phase (SP): in the particular case α = β and both α, β < k/2, both low density and high density regions coexist within the lattice. They are then separated by a domain wall which is non-localised, as it can be shown to diffuse randomly throughout the entire lattice [45].

All these phases are represented in Fig 1a, which constitutes the ‘phase diagram’ of this simplest version of the TASEP. These results can be most easily obtained from the mean-field approximation, but they are known to become exact in the limit of an infinitely long lattice [46].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Main characteristics of the TASEP model.

(a) Phase diagram. A low density (LD) phase, with density ρ < 1/2, arises when the current is limited by the input rate α. Symmetrically, a high density (HD) phase, with ρ > 1/2 is observed when the flow is governed by the output rate β. When neither is limiting, in the ‘maximal current’ phase (MC), the current is as large as the hopping process in the bulk permits (J = k/4 for ρ = 1/2). The line separating the LD and HD phases is where a shock phase (SP) is observed, where LD and HD zones coexist along the lattice. In this case the (few) sites corresponding to the interface between these these two zones is known as the ‘domain wall’. (b)-(e) Density profiles ρ_i of the original TASEP for a lattice of N = 500 sites and hopping rate k = 1. (b): LD regime: α = 0.2, β = 1; (c): HD regime: α = 1, β = 0.1; (d): MC regime: α = β = 0.5; (e): SP regime: α = β = 0.3. The solid black line represents the analytical solution from the mean-field approximation, and the blue points correspond to numerical simulations using the Gillespie algorithm [41].

https://doi.org/10.1371/journal.pcbi.1005555.g001

The corresponding characteristic density profiles (time-averaged occupation numbers along the lattice) are illustrated in Fig 1b–1e. It is thus important to note that the densities given in (i)-(iv) above apply to the bulk of the lattice. Close to the boundaries deviations from these values arise, since the entry and exit rates imply constraints on the boundary sites which must be matched by the density profile. The deviation within these boundary layers can be estimated from the mean-field approach by considering current conservation along the lattice [46].

Incorporation of ribosome drop-off.

As discussed in the introduction, we set out to study the effect of ribosome drop-off along the mRNA, and in particular we wish to predict the effects of ribosome drop-off on both ribosome density profiles and translation rate. We adopt the TASEP model and include an additional detachment rate γ at every site i of the lattice (see illustration in Fig 2). Hence, ribosomes can either hop to the next codon at rate k, given the corresponding site is free, or drop off the lattice at rate γ. We shall show that adding this reaction rate can have a significant impact on the translation rate. Before addressing the question of realistic sequences, we first establish insight into how the drop-off rate influences the traffic dynamics within the model outlined above. We derive analytical expressions for the density profiles in the next section, which we then compare to stochastic simulations of the process, performed using the Gillespie algorithm [41].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Illustration of the ribosome drop-off model.

Ribosomes bind to the first codon of the mRNA with rate α and leave the lattice representing the mRNA strand at the stop codon, with rate β. Throughout the mRNA lattice they can either hop to the next codon, if it is free, with rate k, or drop off the lattice with rate γ.

https://doi.org/10.1371/journal.pcbi.1005555.g002

This first version of the model, based on a uniform drop-off rate, is directly related to the TASEP model with Langmuir kinetics (TASEP-LK), which considers a TASEP of particles which can both bind to and unbind from the lattice at any site [38, 39]. In that case a Langmuir-like equilibrium is established with the environment, and the implications on the current are well understood. However, the discussion in [38, 39] is based on the case where binding dominates with respect to unbinding, a case of interest in the context of transport by molecular motors but not representative of the biology of translation, where only drop-off of particles occurs. Although a mapping between the two cases can be exploited based on mathematical symmetries, it will be important in the following to establish direct intuition for the drop-off problem.

We therefore pinpoint the effects of ribosome drop-off on translation dynamics, which are not obvious in [38, 39]. As a consequence, the analytical expressions derived in our paper are different from the ones in [38, 39], and provide an interpretation for the effects of ribosome drop-off on translation. We therefore adapt our vocabulary to the problem of translation and analyse in detail how different types of mRNAs are affected by ribosome drop-off.

As to the phases, we shall continue to refer to them as LD, HD and MC, even in the presence of drop-off. We do so despite the fact that their characteristics are impacted by drop-off: indeed, since the particle current is no longer conserved along the segment, there is no reason for the density to be constant along the lattice, even within the bulk. In particular, the MC phase no longer corresponds to a density of 1/2 (which, as we will show, is now only true at the left boundary). However, maintaining the nomenclature is a reminder that the nature of the phases remains similar to those of TASEP, in the sense that they are still controlled by the in-rate (initiation rate α) for LD, the out-rate (termination rate β) for HD, and by the bulk hopping rate (elongation rate k) for MC.

Density and current profiles

To derive the density profiles for the stochastic model of ribosome drop-off, we start by writing the master equations in the standard mean-field approximation, i.e., neglecting correlations between neighbouring lattice sites: (5) where we focus the discussion on the steady state, for which .

Eq (5) gives access to the density profile, i.e the local density ρ_i for all sites i along the lattice. However, from a mathematical point of view it is more readily exploited in terms of a continuous function ρ(x) of a positional variable x which varies from 0 to L along the segment. Quite intuitively this amounts to defining the lattice spacing ε = L/N, where N is the number of sites in the lattice. We then require that ρ(x_i) = ρ_i on all lattice points (for which x_i = i ε), and take the ‘continuous’ limit ε → 0 in order to make the lattice points densely cover the interval x ∈ [0, L]. A simple way of stating this is to say that, instead of the discrete lattice, we work with a much more refined one (and in fact an infinitely refined one). We already point out that this implies rescaling some variables [39] in order to make contact with the model predictions, as we will discuss below.

The continuous formulation is thus obtained by identifying ρ_i ≡ ρ(x), ρ_i+1 ≡ ρ(x + ε), ρ_i−1 ≡ ρ(x − ε), and by expanding the density ρ(x) in powers of ε to first order as (6)

For simplicity we will denote ρ(x)≡ρ and from now on, and omit the explicit reference to the spatial dependence.

We first establish the current-density relation J(ρ), by using Eq (6) to express the current J_i = kρ_i(1 − ρ_i+1) in the continuous limit: (7) This shows that, locally, the current still follows relation (4), valid for TASEP, and its variation along the strand is directly due to the density variation. In particular, this implies that the current achieves its maximum value at ρ = 1/2, and we shall see that a density ρ(0) = 1/2 at the entrance of the lattice still defines the MC phase. Moreover, due to relation (7), the LD and HD phases continue to be discriminated by the criterion ρ_LD < 1/2 and ρ_HD > 1/2, respectively. However, in these phases the current is no longer preserved along the lattice.

By substituting Eq (6) into Eq (5) we obtain a condition for the stationary state: (8) Reordering terms yields (9) At this point it becomes clear that, rather than by the individual parameters, the behaviour in the continuous limit is set only by the combination γ/ε. This suggests [39] defining a rescaled drop-off rate as (10) It is this variable which must be matched when comparing predictions to simulations on a discrete system (see also the Supplemental Material for an example). In the context of translation, we are dealing with N codons of a size of about 0.9 nm each, and therefore the rescaled drop-off rate to be used in the continuous model is related to the actual drop-off rate γ by Γ = γ/(0.9 nm).

Consequently, using the rescaled drop-off rate Γ = γ/ε, we thus establish a differential equation for the density profiles as: (11) Notice that Eqs (7) and (11) imply (12) which reflects our intuition about the drop-off process: the current decreases along the lattice due to particles dropping off at each site.

Eq 11 can be solved analytically by separating variables, (13) Notice that the bounds to be applied are specific to the phase and the corresponding boundary conditions. For example, in the LD phase the density at x = 0 is given by ρ_LD(0) = α/k, whereas in the HD phase the density at x = L is ρ_HD(L) = 1 − β/k. Note that in the HD phase it is possible to maintain that value of the density at the right boundary despite the presence of drop-off, since the exit rate is more limiting than the drop-off rate, and therefore a queue of particles can still form. Essentially, all particles which drop off are replenished via the in-rate in this regime. However, there is a critical value of the drop-off rate beyond which the HD phase is not sustainable anymore, as we will see in the section ‘Phase Diagram’. Based on these boundary conditions we can then obtain the density profiles for the LD and the HD phases. Using ρ_LD(0) = α/k, we get (14) which leads to an implicit equation for the density profile: (15) Similarly, by using ρ_HD(L) = 1 − β/k, we get (16) and therefore: (17) Finally, in the maximal current (MC) phase, we expect the current at the left boundary to be close to k/4, the TASEP value, since no significant amount of detachment can have occurred. Therefore, using J_MC(0) = k/4 = α(1 − ρ(0)), the density should be at the left boundary. However, we cannot use ρ(0) as a bound for the integral in Eq (13): the continuous approach, in the above approximation, does not give access to the variation of the density profile within boundary layers. Instead we use the condition J_MC(0) = k/4 = kρ(0)(1 − ρ(0)), which yields ρ_MC(0) = 1/2. Using this bound for the integration, we obtain (18) which leads to (19) Eqs (15), (17) and (19) for the density profiles ρ_LD(x), ρ_HD(x) and ρ_MC(x) have an explicit solution in terms of the so-called Lambert function [38], denoted by W and defined as the inverse function of y = f(z) = ze^z, i.e., z = f⁻¹(y) = W(y). To see this one can multiply both sides of Eqs (15), (17) and (19) by −2, which yields an equation of the form z(x)e^z(x) = y(x), where z(x) = −2ρ(x) and y(x) stands for -2 times the right-hand side of Eqs (15), (17) and (19), respectively. Therefore the solution is given as ρ(x) = −W(y(x))/2.

The W function is multivalued, involving two main branches denoted as W₀ and W₋₁ (see Fig 3). Hence, for the equations above to be well defined we need to identify the appropriate branch of the Lambert function in each case, i.e. for the LD, HD and MC phases. To do this we consider the following: since in the LD and MC phases the particle density on the lattice is between 0 and 0.5, the branch of interest in these two cases is W₀(y) (blue thick line), and specifically the region marked with the solid blue line in Fig 3. Analogously, since the particle density in the HD phase is between 0.5 and 1, the branch of interest is W₋₁(y) (red line), in the region where it highlighted as the solid red line in Fig 3.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Both real branches W₀(y) (bold blue line) and W₋₁(y) (red line) of the Lambert function.

The solid part of the blue line corresponds to the low density and maximal current phases, and the solid part of the red line corresponds to the high density phase.

https://doi.org/10.1371/journal.pcbi.1005555.g003

Hence, an analytical expression can be given for the density profile in the LD regime: (20) where (21)

The HD phase can be discussed analogously. The density profile is given by (22) where (23) And, finally, in the MC phase the density profile is given by (24) where (25) The expressions for ρ_MC, ρ_LD and ρ_HD depend on the ratio Γ/k, as well as on α/k and β/k in the LD and HD cases, respectively. We thus define the rescaled parameters , and and use them instead of the original parameters from now on.

Fig 4 shows three density profiles, alongside the corresponding current profiles, for the LD (a,b), HD (c,d) and MC (e,f) phases. The analytical expressions obtained above (black lines) show good agreement with the numerical simulations (blue points), apart from the boundary layers. Note that this disagreement at the boundaries decreases with increasing number of lattice sites [39].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Density and current profiles numerically computed via stochastic simulations (blue and magenta dots, respectively) and analytically estimated (black lines) for and L = 1: (a, b) LD: = 0.2, = 1; (c, d) HD: = 1, = 0.2; (e, f) MC: = 1, = 1.

Hopping and drop-off events are scheduled based on the Gillespie algorithm [41]. All simulations, here and in the following, are for 10⁷ iteration steps (each iteration corresponds to one reaction: initiation, elongation, drop-off or termination). The first 2 ⋅ 10⁷ iterations were discarded to make sure that the system was in a steady state. Unless otherwise stated, a lattice size of 500 sites was used.

https://doi.org/10.1371/journal.pcbi.1005555.g004

In the LD case, the density profile decreases from left to right (Fig 4(a)). This is as expected: the entry rate fixes the density at the left boundary, and due to the drop-off rate some particles are lost on their journey to the 3’ end (right-hand boundary) of the lattice, causing a decreasing density profile. For the HD phase, however, the ribosome density profile counter-intuitively increases from left to right (5’ to 3’ end) (Fig 4(c)). One could expect that, since more and more ribosomes drop off as they travel along the mRNA, their density should always decrease along the strand. However, in the HD phase there is another factor: as drop-off progressively depletes the ribosome population, their motion is less impeded, which increases the current with respect to the case without ribosome drop-off (see Fig 4(d); the corresponding current without ribosome drop-off would be equal to J(L) and constant throughout the mRNA (not shown)). This constitutes a compensatory mechanism, which ultimately leads to a build-up of particles downstream the segment, i.e. an increasing density profile. It is the balance between these two facts which sets the density slope. Mathematically, this can be seen by considering Eq (12), which implies that J(x) is a non-increasing function of x. Moreover, J(x) is a symmetric function of ρ with a maximum at ρ = 1/2. Therefore, if ρ(0) > 1/2, like in HD, then ρ(x) must be increasing to allow J(x) to decrease. Analogously, if ρ(0) < 1/2, like in LD, ρ(x) must be decreasing (we thank the anonymous reviewer for providing this argument).

A further important characteristic of mRNAs sustaining an HD phase is that, also somewhat counter-intuitively, ribosome drop-off induces an increase in the current at the 5’ end of the mRNA. This is simply due to the fact that the HD phase is limited by the termination rate: therefore current and density at the 3’ end are directly set by the termination rate, at the TASEP values. However, the current decreases along the strand, as discussed above, implying that the current lies above the TASEP value anywhere else in the segment. Summarising, on mRNAs carrying a HD phase, ribosome drop-off induces a current increase at the 5’ end.

Localised domain wall

Much useful insight into the phase diagram can be gained from understanding the circumstances under which a shock phase arises, i.e. when different phases co-exist in upstream and downstream zones. In the presence of drop-off, the position of this domain wall is localised, i.e. the associated density mismatch arises at a specific position along the lattice [39]. This is in contrast to what is seen for the regular TASEP, and we first recall the domain wall phenomenology in this latter case.

In the conventional TASEP, without drop-off, the shock phase (SP) arises between an LD (upstream) zone and a HD (downstream) zone, when α = β. However, the position of the domain wall will fluctuate randomly [45], as is most easily seen by considering the case of very small input/output rates, α = β ≪ 1. In this regime the entry/exit of particles are rare events, and the densities on the respective sections of the lattice are ρ_LD ≃ 0 and ρ_HD ≃ 1 (see ‘Materials and Methods’). We are thus essentially dealing with a queue of particles in front of the exit site, extending up to the domain wall position x_w. Note that this is a valid stationary state, despite the different densities at different sides of the domain wall, since both zones carry the same current (see J_LD and J_HD in ‘Materials and Methods’ for α = β). Consider now an isolated event of an additional particle entering the lattice. Since the exit rate is also very small, the most likely scenario is that the extra particle can reach the domain wall, where it gets stuck. This effectively moves the domain wall upstream by approximately one lattice site, at rate α. Alternatively, the event of a particle leaving the lattice will displace the domain wall downstream by (approximately) one lattice site, at rate β. In both cases the new domain wall position corresponds again to a stationary state, since all currents remain unaltered. Since both domain-wall displacements occur with the same rate (α = β), the domain wall is seen to perform a random walk, and therefore to diffuse freely all over the lattice. This argument can be easily extended to the case of general α = β < 1/2, i.e. it is general to the shock phase in TASEP, and its dynamics are well understood [45].

The presence of drop-off introduces an important difference: the density profiles, and therefore the current profiles, are now position-dependent even throughout the bulk of the LD and HD zones. Nevertheless, the current still has to be continuous across the domain wall: otherwise, the mismatch of currents would locally redistribute particles across the domain wall, which would thereby be displaced until the currents match. Therefore, there is one particular position x_w at which the domain wall can remain, and it is set through the condition J_LD(x_w) = J_HD(x_w), which leads to (26) This is also seen directly from the dynamic picture given above. Consider an entrance event, so that the domain wall would be pushed upstream, say to x_w − ϵ. As is apparent from the current profile, e.g. in Fig 5(f), at this new position the current corresponding to the HD segment is higher than the current of the LD segment. The resulting imbalance therefore leads to the domain wall moving downstream, restoring the initial domain wall position. Using the analogous argument for an exit event, it becomes clear that there is indeed a privileged position at which the domain wall can remain stationary. It is easy to generalise the argument to arbitrary values of α and β, but one may wonder whether the stabilisation mechanism always applies. Indeed, it requires the current slope of the LD phase to be weaker, at the position of the domain wall, than that of the HD phase. But this condition is seen to be always fulfilled by using Eq (12) as well as the fact that ρ_LD(x) < ρ_HD(x) anywhere at the lattice.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Density and current profiles in the shock phase, for , L = 1 and .

(a,b) , (c,d) and (e,f) . Numerical simulations are shown as blue (density) and magenta (current) points, whereas analytical predictions for the density profiles are shown as black lines. For the current profiles the analytical expressions for J_LD(x) and J_HD(x) are shown with green and red dashed lines, respectively.

https://doi.org/10.1371/journal.pcbi.1005555.g005

A closer analysis shows another new feature due to drop-off: there is now a second type of shock, which involves coexistence between the HD and MC phases. Indeed, as we discussed above, in the HD phase the loss of ribosomes along the mRNA leads to a ribosome density profile which increases from 5’ to 3’. In the MC phase, on the other hand, it decreases. Just as in the case of LD-HD coexistence, it is possible to have MC-HD coexistence separated by a domain wall, located at x_w, such that the current remains continuous, i.e. J_MC(x_w) = J_HD(x_w). This condition leads to (27) Importantly, there is a density mismatch at the domain wall position, just as for the LD-HD coexistence phase, and we are thus again dealing with a proper coexistence associated with the MC-HD transition. Note that this contrasts the regular TASEP, where this transition is continuous. Note also that there is no further type of coexistence: using Eqs (20, 21) and (24, 25), together with Eq (7) one can see that the currents in the LD and MC phases coincide only if . At this point the LD and MC phases have identical density profiles, which excludes the scenario of a domain wall and therefore implies a second order transition.

In the following we refer to a phase featuring a domain wall as a ‘shock phase’ (SP), and we only distinguish between a LD-HD and a MC-HD shock phase when this is useful. According to the nature of the interface, either Eq (26) or Eq (27) must hold across the domain wall due to current conservation, from which a unique domain wall position x_w follows. It will be useful in the following to realise that the value of x_w indirectly indicates how the shock phase evolves as one of the rates (initiation, elongation or termination) is varied. Indeed, for a shock to be present the domain wall must lie within the bulk of the lattice (x_w ∈ [0, L]). However, a domain wall position reducing to zero (x_w ≤ 0) indicates that the domain wall has moved to the left boundary, which implies that the HD phase is established all along the lattice. Analogously, x_w ≥ L means that the domain wall has moved to the very right boundary, and therefore the LD phase is established across the entire lattice.

We illustrate this in Fig 5, which shows three density and current profiles corresponding to the SP, for identical conditions except for the exit rate which is progressively increased (0.15 in panels a,b; 0.16 in panels c,d; and 0.17 in panels e,f). We see how the position of the shock evolves towards the right boundary (3’ end) as the value of increases. Eventually, when becomes small enough, the domain wall attains x_w = 0, which indicates that the high density domain occupies the entire lattice beyond this threshold. Predictions from the analytical expressions show good agreement with simulation data, apart from the transition zone from low to high density, close to the position of the domain wall. This discrepancy can be attributed to the fact that finite size effects are present: indeed, Fig 6 shows how the transition becomes sharper as we increase the system size.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Illustration of finite size effects in the shock phase, based on the density profile for , , and L = 1.

The blue dots represent the numerical results for N = 500 sites,whereas red dots are for N = 1000 sites. The solid black line corresponds to the analytical solution.

https://doi.org/10.1371/journal.pcbi.1005555.g006

Phase diagram

Having established a complete catalogue of the phases, we can now construct the phase diagram, and comment in particular on those phases where shocks appear. The full phase diagram must specify, for each point ) in the phase plane, the phase corresponding to these parameters.

The most straightforward observation to make is that the line separating the LD and MC phase remains unchanged with respect to TASEP, i.e. it occurs at . This can be seen from Eqs (20, 21) and (24, 25), having in mind that the density in LD must respect the condition ρ_LD < 1/2. As we have seen, the LD density profile decreases along the sequence in the presence of drop-off, and therefore this condition will first be violated at the left boundary (x = 0) as the rate is progressively increased. We can thus evaluate the boundary of the LD phase through the criterion ρ_LD(0) = 1/2, which yields . Recall that this transition from LD to MC is continuous, as already discussed above: this is seen here directly by substituting into Eq (20).

As a second element to the phase diagram, we recall that the ‘shock phase’ (SP) now englobes two different kinds of coexistence, LD-HD (for ) and MC-HD (for ). They too are therefore separated by a line at .

The remainder of the phase diagram can now conveniently be established by determining the boundaries of the SP, i.e. by determining the conditions at which coexistence ceases. As discussed in the previous section, the SP can be identified as those points for which either Eq (26) or Eq (27) is fulfilled, recalling that x_w must lie within the interval [0, 1]. Equating x_w = 0 amounts to calculating the transition from HD to SP (since the HD has taken over the entire mRNA lattice when the domain wall is located at the left boundary). Analogously, equating x_w = L corresponds to calculating the SP to LD transition. Here we can invoke our insight established above and illustrated in Fig 5: at a fixed we expect the domain wall to be pushed towards the entrance (5’ end) as is reduced, and towards the exit (3’ end) as is increased. Therefore there are both a lower bound and an upper bound to the SP zone. We can determine these by solving for x_w = 0 and x_w = L, respectively, distinguishing the cases (for which the SP coexistence is LD-HD), and (for which it will be LD-MC).

Specifically, the upper bound of the SP can be determined by setting x_w = L in Eq (26) for , and in Eq (27) for , and then solving for as a function of . This yields (28) which defines the LD-SP and MC-SP phase boundaries, respectively. Notice that the transition between both boundaries is continuous at , as expected.

In the same manner, the SP-HD boundary is determined by setting x_w = 0 in Eq (26) for , and in Eq (27) for , and solving for as a function of : (29) which again yields a continuous transition at .

The conditions for the different phases can therefore be summarised as follows:

• LD phase: and .
• MC phase: and .
• HD phase: .
• SP phase:

Alternatively, the maximum current principle [44] can be invoked to deduce the LD-MC and HD-MC boundaries. The maximum current principle states that the value of the current in the MC phase is given by J_MC = max_{[ρN+1, ρ0]} J(ρ), where ρ₀ and ρ_N+1 represent densities associated to particles reservoirs on the left and right hand side of the lattice, respectively. Since in our ribosome drop-off model the fundamental relationship J(ρ) = kρ(1 − ρ) remains valid, the maximum in the current also occurs at ρ = 1/2. Then, by realising that the maximum value of ρ in LD can only be reached at x = 0, and equating ρ_LD(0) = 1/2 we obtain the value of at which the LD-MC transition occurs, namely . Analogously, by realising that the value of ρ that maximises the current in HD can only occur at x = 0, one can equate ρ_HD(0) = 1/2, which yields . Note that these values for where the boundaries occur are consistent with the ones obtained by the previous approach.

In order to validate this analysis we have used simulation data to calculate numerically the average density on the lattice, sampling points ( across the phase diagram. This average density is shown, together with the phase boundaries obtained from the mean-field approach, in Fig 7. As expected, the average density calculated from simulation clearly distinguishes between the LD, HD and MC phases (the latter showing a constant average density across the entire phase, also as expected). Furthermore, to numerically validate the presence of a domain wall we have calculated the value at which the density ρ(x) crosses 0.5 from below (after applying a moving average filter to the simulated density profile spanning 50 codons): this indicates a domain wall, and cannot happen in a LD, HD or MC phase. This method thus provides a straightforward tool to numerically determine the SP region, and furthermore to determine the position of the domain wall. Fig 8 shows those points ( for which a domain wall has been detected from simulation data. The corresponding zone coincides reasonably well with the mean-field prediction for the SP region in the phase diagram, with a slight disagreement, most clearly present for .

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Phase diagram for and L = 1 showing the average density on the lattice (heat map) numerically computed via stochastic simulations, and the borders among the phases (black lines) estimated analytically within the mean-field approximation.

https://doi.org/10.1371/journal.pcbi.1005555.g007

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. SP phase in the phase diagram, for and L = 1.

The colour map represents the value of x_w for the position of the domain wall between LD and HD zones, as determined numerically. Where no data is shown x_w was negative or larger than L, indicating respectively a LD or a HD state point. The represented points thus identify the zone corresponding to a SP, which is to be confronted to the analytical estimation of its phase boundaries (black lines).

https://doi.org/10.1371/journal.pcbi.1005555.g008

To complement the previous analysis, for which a particular value of was chosen, we now show how the extent of the SP region depends on . To this end Fig 9 represents the boundaries of the SP region for various values of , superposed using different colours. The TASEP case without drop-off (), for which the SP zone collapses to a line, is shown as a reference (black line). Initially, as increases, the SP region widens and is progressively shifted to lower values of β (blue and green lines). During this process the LD, MC and SP regions expand at the expense of the HD region. The HD zone ultimately vanishes completely as the lower line delimiting the SP zone is shifted downwards, out of the represented portion of the plane.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Extent of the MC-HD shock phase for L = 1, depending on (analytical results); black line: , blue lines: , green lines: and the purple line: .

The insets show the full phase diagram for these three values.

https://doi.org/10.1371/journal.pcbi.1005555.g009

A complementary view of this progression is shown in Fig 10, where we quantify the extent of the SP zone which involves an MC-HD interface () through its upper bound (blue line, ) and its lower bound (green line, ) as a function of . In the absence of drop-off () both bounds are equal to 0.5, reproducing the behaviour known for TASEP. Then, as increases, the zone between and drops and widens, until reaches zero: this point represents the ribosome drop-off rate at which the HD phase becomes entirely unsustainable on the lattice. Using Eq (29) this value can be calculated as (30) indicated in Fig 10. Once the threshold has been crossed, the SP is delimited by the condition . However, since the upper boundary of the SP asymptotically tends to 0 as , the extent of the SP region progressively vanishes for , with the LD and MC regions taking over the phase diagram: both HD and coexistence phases become inviable if the drop-off rate is too high.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Analytically computed upper and lower boundaries to the MC-HD shock phase, depending on for L = 1.

The green line represents the lower boundary and the blue line corresponds to the upper boundary . The horizontal red dashed line serve as illustration. The intersections (at ) with and indicate the values of the drop-off rates at which an mRNA with termination rate would undergo a transition from HD to SP, and from SP to MC, respectively.

https://doi.org/10.1371/journal.pcbi.1005555.g010

Ribosome drop-off effects depend strongly on codon configuration

The phase diagram is not only key to a systematic characterisation of the model, but is furthermore central to understand biological implications. The link is provided by a genome-wide analysis of ribosome traffic dynamics in S. cerevisiae [17], which has shown that different mRNAs can be classified according to the type of phase transition they present as the initiation rate is increased: in general, mRNAs can either undergo an LD-HD-like transition, or an LD-MC-like transition. Importantly, this classification has also shown a statistically significant correlation with the biological function encoded by the mRNA. For example, mRNAs encoding protein components of the ribosome and translational machinery predominantly undergo an LD-MC-like phase transition as the initiation rate α is increased. In contrast, mRNAs coding for regulatory proteins, such as transcription factors, predominantly undergo an LD-HD-like transition with increasing α. It was found that GO annotations (process) related to stress responses were significantly over-represented in mRNA sequences undergoing an LD-HD-like transition: biological regulation (P-value 1.74 × 10¹¹) and cellular response to stimulus (P-value 1.79 × 10⁻⁸). Moreover, GO annotations (process) related to translation were significantly over-represented in mRNAs undergoing LD-MC-like transitions: cytoplasmic translation (P-value 7.04 × 10⁻¹⁴), translation (P-value 8.58 × 10⁻⁸). The P-value was computed using the hypergeometric distribution, and the Bonferroni correction was applied to take into account multiple hypothesis [54].

This prompts the question whether mRNAs in different phases are affected differently by the ribosome drop-off effect. If they react differently to drop-off, then this could indicate a further regulatory mechanism. In order to quantify the impact of ribosome drop-off, we consider the ribosome current at position x = L of the lattice, since this is the quantity that corresponds to the protein synthesis rate (full-length proteins produced per unit time). In particular, we consider the ratio (31) where J₀(L) denotes the protein translation rate of that same mRNA sequence in the absence of ribosome drop-off. We refer to as the ribosome drop-off resilience of a given mRNA, since characterises the extent to which the tranlation current is maintained despite the drop-off, where all other rates are maintained constant. By using Eqs (20), (22) and (24), together with the expression for the current Eq (7), we can then calculate the drop-off resilience in all LD, HD and MC phases.

Fig 11 shows depending on for mRNAs in three different phases: LD (blue), HD (red) and MC (green). The ribosome drop-off rate has a similar impact on both the LD and the MC regime, showing a striking decrease in the translation efficiency with the drop-off rate . The HD phase, however, is crucially different: if an mRNA strand carries an HD phase, then the ribosome drop-off does not lead to a ribosome current loss at first. This also holds for the SP phase (for which the final section of the sequence is also in HD). This may be understood from the observation that in these phases the current is determined by the exit rate, implying that the density at the last site of the mRNA lattice is , irrespective of the presence or absence of ribosome drop-off in the bulk (see Eqs (22) and (23)). This reflects the fact that, since initiation is frequent, and ribosomes which drop off are effectively replaced, the impact of drop-off is compensated in this regime. Therefore, the protein synthesis rate of those mRNA strands whose dynamics is governed by ribosome queueing, i.e. those carrying HD and SP phases, is not affected by ribosome drop-off. The following section will show that this remains qualitatively correct even for more realistic models of the mRNA strand.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Resilience to drop-off for various choices of initiation and termination rates, computed analytically (solid lines) and numerically (dots) for L = 1.

The blue line represents a mRNA sequence in the low density regime, with . The red line corresponds to a sequence starting in the high density regime with , . The maximal current regime is represented by the green line with and .

https://doi.org/10.1371/journal.pcbi.1005555.g011

Fig 11 also shows, however, that this only holds up to a critical value of the drop-off rate, beyond which the drop-off resilience finally decreases. Indeed, the red curve corresponds to an mRNA which shows perfect resilience at first: its translation rate remains entirely unaffected by a small drop-off rate. Closer inspection shows that the mRNA first carries a HD phase (for very small ), but then undergoes a first transition from HD to SP (left vertical dashed line in Fig 11) as the drop-off rate is increased. This is followed by a second transition from SP to MC (right vertical dashed line in Fig 11), beyond which the resilience starts to decay. The sequence of transitions is easily visualised in Fig 10, where the red dashed line indicates the termination rate associated to this particular mRNA. Where this line crosses and (green and blue lines, respectively) indicates the values of at which the HD-SP and SP-MC transitions occur, corresponding to and in Fig 11, respectively.

To calculate the value we use the fact that the termination rate fulfills for this particular value of . Using Eq (29), this occurs at a ribosome drop-off rate (32)

Analogously, we can find the critical value for the ribosome drop-off rate at which the mRNA subsequently switches from the SP to the MC phase, which occurs when . Using Eq 28 one finds (33) Hence, if exceeds the mRNA can no longer sustain a HD phase, not even close to the exit site. This is the point beyond which the translation rate efficiency starts to decrease with the drop-off rate.

Importantly, it follows from this analysis that mRNAs with a slower termination rate are less susceptible to ribosome drop-off. This is illustrated in Fig 12, which represents the drop-off resilience for three mRNAs with different values of , all of them starting in HD. We indicate with dashed lines the mean-field results and the mean-field critical values and for each of the mRNAs. Simulation results clearly confirm the fact that the translation rate remains unaffected by small drop-off rates, and mRNAs with large exit rates are first to suffer a drop in translation efficiency. As decreases, however, the disagreement between the mean-field expression and the simulations worsens. Increasing the transient and total integration time does not yield any improvement, and increasing the lattice size to N = 1000 leads only to a slight improvement. Simulations show that as decreases, correlations become important, and therefore the mean-field fails to correctly locate the value of at which the SP-MC transition occurs.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 12. Resilience to drop-off, computed analytically (solid lines) and numerically (dots) for L = 1 and .

Red: , blue: , and green: . For these simulations, the transient time was 10⁹ iterations, followed by an integration time of 10⁹ iterations.

https://doi.org/10.1371/journal.pcbi.1005555.g012

In conclusion, our model indicates that the translation rate of mRNAs carrying a HD phase are substantially more resilient to ribosome drop-off compared to other mRNAs carrying LD or MC phases.

mRNAs with high elongation rates are more resilient to ribosome drop-off

We next analyse how the resilience of mRNAs to ribosome drop-off is affected by the elongation rates k. We start by computing the resilience of a homogeneous mRNA sequence depending on its elongation rate k. We have seen in the previous section that mRNA sequences can be classified according to the phase they can carry, and the translation rate of mRNAs that sustain an HD phase is not affected at all by ribosome drop-off, i.e., they achieve the maximal value of resilience. We therefore analyse the effect of the elongation rate on the resilience for 3 different mRNA sequences, each one carrying an LD, MC or HD phase, respectively. The results are shown in Fig 13, where the solid lines show the analytical results and the dots represent the results from the numerical simulations. It is clearly seen that for both LD and MC, the resilience to drop-off increases monotonously with the elongation rate k. In contrast, for the sequence carrying HD, the resilience remains constant at the maximal possible value, χ_Γ = 1, as expected. Notice that the range of values of the elongation rate k for each curve is chosen so that the corresponding phase (LD, MC or HD) is sustainable with the given parameters. For example, for the chosen parameters shown in Fig 13, to be in LD we have to fulfil α < k/2; the condition for MC is α > k/2, and for HD we must ensure that β < β_l, with β_l defined as in Eq (29) (where one must take into account that β_l depends on k).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 13. Drop-off resilience versus elongation rate analytically (solid lines) and numerically computed (dots) for L = 1 and Γ = 0.05.

Blue (LD): α = 0.1, β = 1; Red (MC): α = 1, β = 1; Black (HD): α = 1, β = 0.1. For these simulations, the transient time was 10⁷ Gillespie iterations, followed by an integration time of 10⁸ iterations.

https://doi.org/10.1371/journal.pcbi.1005555.g013

The results obtained for the LD and MC phases are quite intuitive: high values of elongation rate lead to fast ribosome transit through the mRNA, and consequently, the probability that they are affected by drop-off is smaller. Hence, higher elongation rates lead to higher resilience. However, the results for the HD phase might appear counter-intuitive; since in HD ribosomes have longer dwell times due to queue formation, one might expect lower values for the resilience compared to the LD and MC cases. This can be understood, however, realising that in HD the translation rate is determined by the ribosome density at right boundary of the lattice. Due to the fact that the exit rate in this phase is more “limiting” than the drop-off rate, ribosome drop-off does not affect the translation rate value. Nevertheless, the probability for a ribosome to finish translation is indeed smaller in the HD phase than in the other phases (see ‘Supplementary Information’).

The numerical results for the MC phase show some discrepancy with respect to the analytical results. This is due to the higher correlations present in the MC phase with respect to the other phases, neglected in the mean-field approximation. Moreover, there is also a slight discrepancy between analytical and numerical results in the minimal value of the hopping rate k needed to sustain HD for the given α, β, and Γ parameters (β < β_l, (Eq (29)). This indicates that there is a slight discrepancy between the mean-field prediction for the boundaries between HD and SP, as seen in Fig 8.

Taken together, these results suggest that mRNA sequences using fast codons are more resilient to ribosome drop-off than mRNA sequences containing slower codons. This is consistent with the long-standing observation that mRNAs whose usage of different synonymous codons is highly biased (using predominantly fast codons) code for proteins needed at a high level, such as ribosomal proteins and proteins involved in translation [47–49]. Hence, it suggests a regulatory mechanism by which such mRNAs not only have a higher translation rate, but they are also less vulnerable to ribosome drop-off.

Application to real mRNA sequences

We now apply the ribosome drop-off model to mRNA sequences from S. cerevisiae. As mentioned previously, for such real mRNAs with a specific codon sequence, the elongation rate is in fact codon dependent. This amounts to inhomogeneous hopping rates along the lattice, according to the mRNA sequence. Here we analyse to which extent the conclusions drawn from the theoretical analysis in the previous sections, essentially based on mono-codon sequences, hold for real mRNA sequences.

To better match our model to the context of translation, we also account for the size of the ribosomes, by acknowledging that a ribosome should be taken to cover 9 codons [7]. In the language of our lattice model, this means that a particle (ribosome) now covers 9 lattice sites, although it still steps ahead by single sites. This question has been explored in several studies [50, 51]. Although it leads to a significant modification in the current-density relation, the qualitative phenomenology remains very similar. This refinement therefore serves to better represent the physiological conditions, but should not compromise the conclusions established above.

We have determined the hopping rates k_i associated to each of the codons in the sequence based on the corresponding tRNA abundance for S. cerevisiae, and used the estimated termination rate β = 18.03 s⁻¹, as detailed in ‘Supplementary Information’. The value of the termination rate β is assumed to be the same for all mRNAs, since it depends on the concentration of available release factors. We now analyse the effect of two different values of the ribosome drop-off rate estimated in [29]: γ = 1.4 × 10⁻³ s⁻¹, i.e. the value that has been estimated under physiological conditions, and γ = 5.6 × 10⁻³ s⁻¹, estimated for ethanol intoxication conditions, both for E. coli. We then compare with results without drop-off (γ = 0) for reference. In [29] the values for the drop-off rate are given in units of ‘events/codon’. We transform them to units of ‘events/second’ by considering that the average speed of ribosomes is 10 codons/s. For example, 1.4 × 10⁻⁴ events/codon means that on average one ribosome needs to travel 1/1.4 × 10⁻⁴ codons before it drops off the mRNA. Equating this distance to the average ribosome speed times the average time between drop-off events allows calculating the latter, and therefore we obtain the frequency as the inverse of the average time between drop-off events. To our knowledge, there are no measurements of ribosome drop-off rate for yeast, and we thus adopt the same rates as for E. coli, since the mechanism of translation elongation is similar in prokaryotes and eukaryotes, and E. coli and yeast are both rapidly growing microorganisms found in environments where competition demands similar optimisation of the translation system [24, 25, 27, 28]. Moreover, we point out that we do not expect our results to change qualitatively even in the case that the drop-off rate γ turned out to be significantly different from the one estimated for E. coli, based on the results presented in Fig 11: as the resilience decreases monotonously with γ, all phases would be affected in a similar way. The only scenario for a qualitative change to occur would be for mRNAs containing an elongation bottleneck. For these, a substantially higher value of γ might preclude them from sustaining a HD phase, and it is only for mRNAs exhibiting a very severe elongation bottleneck that a HD regime might survive.

Our theoretical results in the previous section indicate that the impact of these ribosome drop-off rates strongly depends on the phase of the mRNA sequence under consideration: low density, high density or maximal current regime. But additionally, the effects of ribosome drop-off will clearly be stronger the longer the mRNA sequence is. The rescaled drop-off rate Γ = Nγ/L served to take this into account, precisely by making the results comparable whatever the length of the sequence. Here however, we focus on the effect of ribosome drop-off on translation in real mRNAs, and we thus reason in terms of the actual (non-rescaled) drop-off rate γ. We thus select three representative mRNA sequences from S. cerevisiae, all of them consisting of approximately 500 codons, and corresponding to the different regimes LD, HD and MC. For each sequence we choose the appropriate value of the initiation rate to ensure that it is in the corresponding phase (LD, HD and MC, respectively). These values do not correspond to the physiological initiation rate values; this simulation just illustrates the effects of drop-off depending on the phase for real mRNAs. In the next section we analyse the effect of drop-off for the S. cerevisiae genome under physiological values of the initiation rate.

We first analyse the mRNA sequence of CDC7 (YDL017W), involved in origin firing and replication fork progression. At the initiation rate value of α = 7 × 10⁻² s⁻¹ this mRNA is expected to carry an LD phase. Fig 14 shows the results for the density and current profiles for this sequence, which clearly qualify as LD-like, subject to the expected disorder along the sequence. As the ribosome drop-off rate increases (middle and lower panels) we obtain density and current profiles which decrease along the lattice, as predicted by the theoretical results in the section ‘Density and Current Profiles’, and indeed the slope grows more important for higher γ. Moreover, the current at the first site of the mRNA remains constant as we increase γ, also in accordance with the theoretical results.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 14. Density (a,c,e) and current (b,d,f) profiles for CDC7 (YDL017W), for α = 7 × 10⁻² s⁻¹, in a low density-like phase.

In the density profiles we represent the coverage density of ribosomes, i.e., the probability for a codon to be covered by a ribosome. (a),(b): no ribosome drop-off γ = 0; (c), (d): physiological value of ribosome drop-off γ = 1.4 × 10⁻³ s⁻¹; (e), (f): ribosome drop-off under ethanol stress γ = 5.6 × 10⁻³ s⁻¹. The slope in the current profile clearly increases with the ribosome drop-off rate, and in the density profile the slight slope towards the 3’ is also reinforced by drop-off.

https://doi.org/10.1371/journal.pcbi.1005555.g014

As a second example we consider the mRNA sequence of SLT2 (YHR030C), a serine/threonine MAP kinase involved in regulating the maintenance of cell wall integrity, cell cycle progression, and nuclear mRNA retention in a heat shock. For a value of α = 1s⁻¹ this sequence is in a high density-like regime. The corresponding density and current profiles are shown in Fig 15. The ribosome queue is clearly visible in the upper panels with no ribosome drop-off (γ = 0), where the density is almost 1 along the queue, before it drops dramatically after approximately codon 400, where the bottleneck in translation is caused by a cluster of slow codons. Then, for the physiological value of γ (Fig 15c), the density profile in the queue shows a slightly steeper slope from the 5’ to the 3’ end, reminiscent of the increasing density profiles obtained for HD sequences in the section ‘Density and Current Profiles’. The ribosome queue then starts to vanish for higher values of γ, corresponding to ethanol stress, as illustrated in Fig 15(e). Importantly, the current at the end of the lattice remains constant as we increase γ (compare Fig 15b and 15d), as predicted for sequences in HD, so that the current at the beginning of the mRNA increases with respect to the case without ribosome drop-off (compare panels b, d and f at the start codon, i = 0). This implies that the increased drop-off rate leads to a higher current at the beginning of the mRNA, again characteristic of mRNAs in HD, as shown in the section ‘Density and Current Profiles’. Furthermore, we see that for this sequence, the critical value of the drop-off rate beyond which the resilience starts decreasing is below the physiological γ, even though there are still clear signatures of an HD-like phase for the physiological γ. We expect this critical value of γ to be mRNA-dependent, determined by the severity of the elongation bottleneck present in the mRNA.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 15. Density (a, c, e) and current (b, d, f) profiles for SLT2 (YHR030C), for α = 1s⁻¹, a in high density-like phase.

Details, including the values of γ, are as in Fig 14. (a),(b): no ribosome drop-off γ = 0; (c), (d): physiological value of ribosome drop-off γ = 1.4 × 10⁻³ s⁻¹; (e), (f): ribosome drop-off under ethanol stress γ = 5.6 × 10⁻³ s⁻¹. The increasing slope in the 5’ to the 3’ direction characteristic of the HD phase is clearly visible in panel (c). Moreover, panel (d) illustrates the marked increase in the current at the 5’ end, compared to the case without ribosome drop-off (see panel (b)), also characteristic of the HD phase.

https://doi.org/10.1371/journal.pcbi.1005555.g015

Finally, we analyse the mRNA sequence of CRH1 (YGR189C), a chitin transglycosylase involved in chitin transfer in the cell wall. For a value of α = 1s⁻¹ this sequence is in a MC-like regime. Fig 16 shows the effects of increasing ribosome drop-off on the density and current profiles. The current profiles (Fig 16b, 16d and 16f) clearly show a steeper slope as γ increases, qualitatively similar to the predictions for sequences in MC in the section ‘Density and Current Profiles’. The effects on the density profiles (Fig 16a, 16c and 16e) are not as obvious as on the current profiles, due to the inhomogeneous hopping rates. However, a decrease in the density can be observed when focussing on the 3’ end of the mRNA.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 16. Density (a, c, e) and current (b, d, f) profiles for CRH1 (YGR189C), for α = 1s⁻¹, in a maximal current-like phase.

Details, including the values of γ, are as in Fig 14. (a),(b): no ribosome drop-off γ = 0; (c), (d): physiological value of ribosome drop-off γ = 1.4 × 10⁻³ s⁻¹; (e), (f): ribosome drop-off under ethanol stress γ = 5.6 × 10⁻³ s⁻¹. The current profiles clearly show an increase in the slope as ribosome drop-off increases, qualitatively similar to the effects predicted for sequences in MC.

https://doi.org/10.1371/journal.pcbi.1005555.g016

In order to quantify the effect of ribosome drop-off on each sequence, we calculate their individual drop-off resilience (defined in Eq 31), for each of the considered values of γ. The results are summarised in Table 1. As predicted by the theoretical results, we see that the mRNA in the HD-like phase, SLT2, is the least affected by ribosome drop-off: the loss is only 5% in translation rate under physiological conditions, compared to 14% and 13% for the other two mRNAs in LD and MC-like phases, respectively. The advantage of the HD-like phase with respect to ribosome drop-off is also clear under stress conditions, with a loss of only 27% in translation rate compared to 46% and 43% for the other two phases.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Table summarising the results for the drop-off resilience (defined in Eq (31) for the 3 mRNA sequences CDC7 (LD), SLT2 (HD), and CRH1 (MC), and for three different values of γ.

https://doi.org/10.1371/journal.pcbi.1005555.t001

Therefore, these results strongly suggest that the mechanism which has emerged from our theoretical analysis based on modelling single-codon mRNA is robust: also for real mRNA, the impact of drop-off rates is seen to be less important for HD-like sequences.

Genome-wide analysis

Applying the ribosome drop-off model to 3 real mRNA sequences has illustrated the effect of drop-off, showing that its impact depends on the ribosome traffic regime (LD, HD and MC). We now evaluate the resilience to drop-off for the S. cerevisiae genome under physiological conditions and analyse how the biological function encoded by the mRNAs might be related to its resilience.

We start by estimating the physiological value of the initiation rate for each mRNA sequence by using genome-wide experimental data from [52] of ribosome density ρ_ϕ under physiological conditions (see ‘Supplementary Information’). For each mRNA sequence we then identify the physiological value α_ϕ of the initiation rate, using our ribosome drop-off model by requiring ρ(α_ϕ) = ρ_ϕ, i.e. we take the initiation rate to be the one which replicates the experimentally measured ribosome density (following the approach in [17]). The histogram of the α_ϕ values obtained in this way is shown in the ‘Supplementary Information’. With those values of α_ϕ it is possible to compare the model predictions for the translation rate to experimental data. To do this, we have calculated the Spearman’s rank correlation coefficient between and P, where M and P denote the experimentally measured mRNA and protein levels, respectively, taken from [53], and denotes the translation rate estimated with our model. We obtain a value of 0.57, p-value< 10⁻¹⁴, yielding an improvement compared to other protein level predictors, such as the tAIc (Spearman’s rank = 0.38, p-value<10⁻⁶) [17, 53].

Using the ribosome drop-off rate γ = 1.4 × 10⁻³ s⁻¹ as estimated under non-stress conditions [29], we can then exploit Eq (31) to compute the ribosome drop-off resilience for each mRNA sequence. The histogram of the resilience values χ_γ is represented in Fig 17. We see that the effect of the drop-off rate on the translation depends very much on the specific mRNA sequence, with values of the resilience ranging from 0.2 to approximately 1. Therefore, ribosome drop-off can lead up to an 80% decrease in the translation rate under physiological conditions, depending on the parameters characterising the mRNA sequence. Nevertheless, the mean resilience is 0.87 and the median resilience is 0.89, indicating that most mRNAs are only moderately affected by ribosome drop-off under physiological conditions.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 17. Histogram of ribosome drop-off resilience values χ_γ for the S. cerevisiae genome under physiological conditions.

The colours of the different regions correspond to the following χ_γ intervals: green: χ_γ < 0.7; yellow: 0.7 ≤ χ_γ < 0.8; orange 0.8 ≤ χ_γ < 0.9, and red: χ_γ ≥ 0.9.

https://doi.org/10.1371/journal.pcbi.1005555.g017

We next investigate whether there is a correlation between the resilience to ribosome drop-off of an mRNA and the biological function of the protein it codes for. To address this question we divide the histogram into four regions, from small to high values of χ_γ, and we run a GO enrichment analysis on each region (summarised in Table 2). This analysis was done using the software Gene Ontology Term Finder (http://go.princeton.edu/cgi-bin/GOTermFinder) developed at the Lewis-Singer Institute at Princeton. P-values are computed using the hypergeometric distribution, and the Bonferroni correction is applied to take into account multiple hypothesis tests [54]. We thus identify, for each of these regions, the GO categories which are significantly over-represented. The following observations emerge:

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 2. Table summarising the results for Gene Ontology (GO) enrichment analysis for the 4 different regions of the histogram of the ribosome drop-off resilience values (Fig 17) obtained for the S. cerevisiae genome.

https://doi.org/10.1371/journal.pcbi.1005555.t002

In the two regions which correspond to the mRNAs with the lowest ribosome drop-off resilience (coloured green and yellow, respectively), the dominant GO categories are related to biological regulation, response to stimulus, regulation of metabolic process, signal transduction and cell cycle. In the next region (orange), with a fairly high value of drop-off resilience, the salient GO categories are related to central metabolism processes. Finally, in the region with the strongest resilience against drop-off (red), a highly disproportionate number of mRNAs code for proteins involved in cytoplasmic translation, ribosome assembly, as well as rRNA export and transport. Hence, in essence, mRNA sequences coding for translation and metabolism-related proteins are predominant amongst those sequences with high resilience, whereas mRNAs with low drop-off resilience mainly code for regulatory proteins and those involved in stress responses. A GO enrichment analysis taking the same number of genes from each of the four highlighted histogram regions yields results in agreement with the results obtained with the whole sections (see Supplementary Tables).

It is interesting to link this observation to the fact that mRNA sequences relevant to translation and central metabolism are known to be abundant under physiological conditions, whereas those concerning stress-regulation are present at much lower levels only. This resullt is consistent with the finding in [22] that the ribosome drop-off effect is weaker for mRNAs coding for highly expressed genes, even though their analysis was based on only 10 genes. There, the authors argue that highly expressed genes consume more ribosomes, due to their higher mRNA levels, and that it is therefore beneficial for the cell to optimise their codon sequence in order to minimise losses in translation due to ribosome drop-off. This analysis therefore raises the question whether resilience is an important feature of highly expressed genes, which would point to an evolutionary selection mechanism operating to optimise their resilience to ribosome drop-off.

Identifying the biological functions which the classes of mRNA predominantly refer to also provides an indication of the origin of the high resilience of translation-related mRNA sequences. Indeed, these mRNAs have been shown to be highly codon-biased, i.e. they have been optimised to ensure high elongation rates [47–49]. From our mathematical analysis in the previous section, we know that high elongation rates help to minimise the effect of drop-off (see Fig 13). The argument thus points to the bias towards fast codons as the origin for high resilience.

All in all, the analysis on a genome-wide scale reveals several properties which appear to be linked to resilience and attributes a coherent set of properties to mRNAs with specific gene ontologies. One one hand, as already stated, the value of the resilience χ_γ correlates with the elongation rate k, in agreement with our observation above that high elongation rates help to reduce the impact of drop-off. On the other hand, the resilience χ_γ also correlates with the physiological value of the initiation rate α_ϕ (Pearson correlation coefficient R = 0.42, p-value p = 3.36 × 10⁻²³²), and the physiological value of the translation rate (R = 0.58, p-value p < 2.2251 × 10⁻³⁰⁸). Interestingly, for sequences in LD phase, the resilience to drop-off slightly decreases as α increases (see Supplementary Information), but the positive effect of high elongation rates on the resilience χ_γ is stronger. Hence, even though mRNAs using fast codons tend to have high initiation rates, they exhibit a higher resilience to ribosome drop-off compared to mRNAs containing slow codons.

All these observations support the fact that proteins involved in translation and central metabolism are known to be highly expressed under physiological conditions, in contrast to regulatory and stress-related proteins. Consequently, mRNAs coding for highly abundant proteins are typically characterised by high elongation and initiation rates, leading to a high translation rate. This is indeed confirmed on the genome-wide analysis of physiological α-rates, accompanied by a GO enrichment analysis (shown in the ‘Supplementary Information’).

From an additional angle it is also interesting to investigate how the length of the mRNA sequences is linked to their resilience. The resilience values χ_γ are highly anti-correlated with the mRNA length (R = −0.97 and p-value P < 2.2251 × 10⁻³⁰⁸). In view of the drop-off model, this behaviour can again be rationalised: the ribosome drop-off mechanisms, and therefore the drop-off rate γ, are essentially the same for all mRNAs. Therefore the impact of drop-off will be larger for longer mRNAs (see more detailed discussion in ‘Supplementary Information’). Accounting for ribosome drop-off might therefore offer an explanation for the long-standing observation that the average ribosome density on mRNAs in anti-correlated with the mRNA length, as experimentally measured in [7, 55], and still under debate [11, 17, 29].