# Physical modeling of ribosomes along messenger RNA: estimating kinetic rates from ribosome profiling experiments with a ballistic model

Carole Chevalier, Jerome Dorignac, Yahaya Ibrahim, Armelle Choquet, Alexandre David, Julie Ripoll, Eric Rivals, Frederic Geniet, Nils-Ole Walliser, John Palmeri, et al.

► To cite this version:

Carole Chevalier, Jerome Dorignac, Yahaya Ibrahim, Armelle Choquet, Alexandre David, et al. Physical modeling of ribosomes along messenger RNA: estimating kinetic rates from ribosome profiling experiments with a ballistic model. PLoS Computational Biology, 2023, 19 (10), pp.e1011522. 10.1371/journal.pcbi.1011522. hal-03832368

HAL Id: hal-03832368

https://hal.science/hal-03832368

Submitted on 27 Oct 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
RESEARCH ARTICLE

# Physical modeling of ribosomes along messenger RNA: Estimating kinetic parameters from ribosome profiling experiments using a ballistic model

Carole Chevalier¹, Jérôme Dorignac¹, Yahaya Ibrahim¹,², Armelle Choquet³, Alexandre David³, Julie Ripoll⁴, Eric Rivals⁴, Frédéric Geniet¹, Nils-Ole Walliser¹, John Palmeri¹, Andrea Parmeggiani¹, Jean-Charles Walter¹*

1 Laboratoire Charles Coulomb (L2C), Univ. Montpellier, CNRS, Montpellier, France, 2 Department of Physics, Faculty of Natural and Applied Sciences, Umaru Musa Yar'adua University, Katsina, Nigeria, 3 Institut de Génétique Fonctionelle (IGF), CNRS, Montpellier University, Montpellier, France, 4 Laboratoire d'Informatique, de Robotique et de Microélectronique de Montpellier (LIRMM), CNRS, Montpellier University, Montpellier, France

* jean-charles.walter@umontpellier.fr

<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

OPEN ACCESS

**Citation:** Chevalier C, Dorignac J, Ibrahim Y, Choquet A, David A, Ripoll J, et al. (2023) Physical modeling of ribosomes along messenger RNA: Estimating kinetic parameters from ribosome profiling experiments using a ballistic model. PLoS Comput Biol 19(10): e1011522. https://doi.org/10.1371/journal.pcbi.1011522

**Editor:** Alexandre V. Morozov, Rutgers University, UNITED STATES

**Received:** March 31, 2023

**Accepted:** September 17, 2023

**Published:** October 20, 2023

**Copyright:** © 2023 Chevalier et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** The programs performing the fit of the Ribo-seq profiles (for homogeneous and inhomogeneous elongation rates) are freely available at https://github.com/CaroleChe/Ballistic_model_mRNA_fit.git.

**Funding:** This project received partial financial support from the LabEx NUMEV (ANR-10-LABX-0020) within the I-SITE MUSE of the University of Montpellier via the Flagship Project "Gene Expression Modeling" (NUMEV Etendard-2016-09).

## Abstract

Gene expression is the synthesis of proteins from the information encoded on DNA. One of the two main steps of gene expression is the translation of messenger RNA (mRNA) into polypeptide sequences of amino acids. Here, by taking into account mRNA degradation, we model the motion of ribosomes along mRNA with a ballistic model where particles advance along a filament without excluded volume interactions. Unidirectional models of transport have previously been used to fit the average density of ribosomes obtained by the experimental ribo-sequencing (Ribo-seq) technique in order to obtain the kinetic rates. The degradation rate is not, however, accounted for and experimental data from different experiments are needed to have enough parameters for the fit. Here, we propose an entirely novel experimental setup and theoretical framework consisting in splitting the mRNAs into categories depending on the number of ribosomes from one to four. We solve analytically the ballistic model for a fixed number of ribosomes per mRNA, study the different regimes of degradation, and propose a criterion for the quality of the inverse fit. The proposed method provides a high sensitivity to the mRNA degradation rate. The additional equations coming from using the monosome (single ribosome) and polysome (arbitrary number) ribo-seq profiles enable us to determine all the kinetic rates in terms of the experimentally accessible mRNA degradation rate.

## Author summary

Genetic translation, involving the transport of ribosomes, is one of the main steps of genetic expression. The molecular mechanisms that regulate ribosome motion along messenger RNA (mRNA) remain to be elucidated. Ribo-sequencing experiments are now
CC acknowledges doctoral thesis support from the French CNRS. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

**Competing interests:** The authors have declared that no competing interests exist.

capable of providing access to the gene-dependent density of ribosomes along mRNA. The physical modeling of ribosome motion is carried out using 1D transport models, such as TASEP-like ones, which enables one to obtain the kinetic (initiation and elongation) rates from the Ribo-seq densities. Currently used methods, however, do not allow access to mRNA degradation rates and previous experiments lack input to perform a self-consistent quantitative analysis. Here, we propose a new experimental and modeling protocol to obtain the kinetic rates (including mRNA degradation) with better accuracy. Ribo-seq is performed after splitting the mRNA population into groups that depend on the number of translating ribosomes (from 1 to 4). We thus obtain four Ribo-seq profiles per gene properly normalized without requiring independent experiments. We show that these low ribosome number profiles display a higher sensitivity to mRNA degradation than the full polysome profile and discuss the biological translation regimes of several organisms.

## Introduction

The translation of messenger RNA (mRNA) is, together with DNA transcription, one of the two main steps of gene expression. The information encoded into mRNAs is translated into a polypeptidic chain through the action of ribosomes processing along mRNA. Ribosomes match a triplet of nucleic acids (a codon) to the corresponding amino acid. Different kinetic rates are involved during translation: the initiation rate, i.e. the rate at which a ribosome assembles at the start codon of mRNA and starts the translation, the termination rate, i.e. the rate at which the ribosome releases the protein and disassembles from mRNA when it achieves translation, and the hopping rates from codon to codon (elongation rates) which depend on the translated codon (and thus on the position along the transcript) [1].

Despite its importance in gene expression, and more generally in the metabolism of the cell, much remains to be learned about the kinetics of translation. The onset of the high-throughput technique known as ribosome-sequencing (Ribo-seq) gives access nowadays to genome-wide ribosome density profiles along mRNA [2]. Ribo-seq is an experimental technique based on the sequencing of ribosome-protected fragments of mRNA after fixation of ribosomes and degradation of unprotected mRNA. Although this approach has revolutionized the study of translation, one still needs a deeper physical understanding of the molecular processes in order to turn Ribo-seq data into quantitative modeling tools. Ribo-seq data could then be used to probe the microscopic parameters tuning gene expression at the translational level and to understand the role of deregulation of translation in diseases like cancer [3].

Several approaches have been proposed to describe quantitatively Ribo-seq data by means of physical modeling involving one dimensional transport [2, 4–10]. However, fitting the kinetic model parameters from Ribo-seq data currently faces several major limitations. Firstly, estimating kinetic rates from ribosome profiles is hindered by the lack of sufficient independent experimental data. Since this is a typical inverse problem [11], parameters can be under/over-determined creating technical difficulties for making reliable fits. Secondly, most quantitative models assume that mRNAs have infinite lifetimes. Thirdly, the normalization of Ribo-seq data, necessary for a physical modeling, requires the use of data from other experiments, thus introducing an extra uncertainty. The question therefore arises concerning the possible influence of finite mRNA lifetime on the evaluation of ribosome kinetic rates [4–6].

In this article, we model the motion of ribosomes along mRNA using a ballistic model first introduced by Valleriani et al. [4] of unidirectional transport along a filament including finite mRNA lifetimes. Our model does not account for excluded volume interactions, thus traffic
jam effects induced in the Totally Asymmetric Simple Exclusion Process (TASEP) [12], the paradigmatic model for translation modeling, are not possible. The ballistic model, however, remains a good approximation to TASEP for long enough mRNA lifetimes in the low density regime, which is the relevant regime for translation [13]. We mention also the Ribosome Flow Model [7], which displays a mean field approach similar to ours. We go beyond previous approaches, however, by modeling a new experimental procedure where the mRNA population is sorted out into *k*-somes having a fixed number *k* of ribosomes, going from 1 to 4, present at the instant of the sequencing. Our goal is to demonstrate how to exploit this experimental data to establish a reliable method for an absolute evaluation of kinetic parameters (i.e., in time units).

Heyer et al. [14] previously presented an interesting pioneering Ribo-seq study on *Saccharomyces cerevisiae* that was restricted to monosomes and polysomes. According to these authors, mRNA translation could in general be divided into three classes. Although they discuss the influence of mRNA lifetime and conclude that their data reveal a relationship between it and the monosome occupancy, they do not propose a mechanistic explanation, but simply note that genes enriched in monosomes had a lower median mRNA half-life than those without enrichment. Despite this finding, they do not attempt to integrate mRNA lifetime into their classification. We present here a more quantitative analysis that reveals the crucial importance of mRNA lifetime in classifying genes. Our main result is that although mRNA lifetime does not affect the polysome footprint, it does have a strong impact on the monosome one, despite the low probability of having a monosome.

The paper is organized as follows. In the first section, we define the ballistic model coupled to an mRNA degradation process. In the second section, we solve it analytically by calculating the exact *k*-some probability distribution and density profiles as functions of the kinetic parameters. In the third section, we perform a parametric analysis of these quantities showing that: 1) the smaller the value of *k*, the more these quantities are sensitive to mRNA degradation; 2) for typical biological values of the parameters, low order *k*-somes are more sensitive to mRNA degradation than the full population of mRNAs (polysomes). In the fourth section, we use the model to analyze real Ribo-seq data obtained from a human histonic gene. For this particular gene, we give—for the first time to our knowledge—evidence of degradation on Ribo-seq data. These results validate our approach and provide motivation for setting up and testing a quantitative method in the fifth section. That method rests on combining *k*-some (monosome actually) data with polysome data while taking advantage of the sharply increased sensitivity of mRNA degradation on low order *k*-somes with respect to polysomes. As a proof-of-concept, we first apply our fitting method to data generated from the ballistic model itself. By using the experimentally accessible degradation rate to fix the timescale, our approach leads to an improved method for obtaining absolute kinetic rates of ribosomes. In the sixth section, we rationalize our findings by proposing a global picture of the different dynamical regimes in the form of a *phase diagram* based on the ballistic model. It presents a condensed global view of the different mRNA degradation dependent regimes and allows at a glance a general qualitative analysis for any given gene.

## Definition of the ballistic model and *k*-somes

We model the translation of mRNA by ribosomes, and the degradation of mRNA by means of a ballistic model with degradation as sketched in Fig 1. The mRNA is modeled by a filament of length *L* where ribosomes are treated as point-like particles i.e. without excluded volume interactions. This is motivated by the low density of ribosomes on mRNA. For biological parameters, the discrepancy between TASEP and the ballistic model is less than 5% in homogeneous
<figure><img src="image_4.png" /><figcaption>Fig 1. Graphical illustration of ribosome transport with the ballistic model including degradation. The mRNA is represented as a red filament of length L. The ribosomes are point-like particles processing along the filament with a position dependent rate p(x) (where x ∈ [0; L] is the position along the filament). Bound and unbound ribosomes are respectively represented by dark blue and light blue discs. A ribosome can bind at the entry of the filament at a rate α. A ribosome located at the exit unbinds at a rate p(L) = β. The mRNAs get degraded at a rate ω.</figcaption></figure>

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

stretches of elongation (See S1 Fig). Note, however, that this discrepancy can be higher in regions where correlations between codons are higher in the case of a combination of a fast codon followed by a slow codon. This situation creates a traffic jam in the TASEP, which is absent in our model. We leave the analysis of these marginal corrections to a forthcoming paper [13].

Particles get attached at the entrance of the filament with a stochastic rate α (initiation rate at the start codon) and subsequently progress in one direction along the filament with a local deterministic velocity p(x) (related to the local elongation rate), that depends on the position x along the transcript. We shall simply use p(x) = p for a constant (codon independent) hopping rate. Particles leave the filament with velocity p(L) (termination rate at the stop codon).

Let $\mathcal{T}(x)$ be the time for a ribosome to go from genomic coordinate 0 to x. Since the ballistic trajectory of a ribosome, $X_0(t)$, obeys $dX_0(t)/dt = p(x)$, we obtain

$$ \mathcal{T}(x) = \int_{0}^{x} \frac{dy}{p(y)}. \qquad (1) $$

The particular value $\mathcal{T}(L)$ corresponds to the time needed for a ribosome to cross the entire mRNA. This quantity plays an important role in the determination of the average density profiles of ribosomes along mRNA.

Our model includes the degradation of mRNA at a constant rate ω as in Ref. [4]. For the sake of simplicity, we assume instantaneous mRNA degradation—translation is therefore instantaneously aborted—but other, more progressive, mechanisms exist [4, 6]. We assume that mRNAs are synthesized and degraded in the cytosol in such a way that their population is constant on average. In what follows, we shall refer to the “age”, a, of an mRNA as the time elapsed between its synthesis and the time at which the Ribo-seq experiment is actually done. The mRNA lifetime θ is modeled by a random variable with probability density function $\phi_{\theta}(\theta)$. The assumption of stationarity of the mRNA population leads to a relation between the mRNA age distribution, $\phi_a(a)$, and the lifetime distribution $\phi_{\theta}(\theta)$ [15]

$$ \phi_a(a) = \frac{1}{\langle \theta \rangle} \int_a^{\infty} \phi_{\theta}(\theta) d\theta, \qquad (2) $$

where $\langle \theta \rangle = \int_0^{\infty} \theta \phi_{\theta}(\theta) d\theta$ is the average lifetime of an mRNA. As in [4], we shall consider an exponential mRNA lifetime distribution $\phi_{\theta}(\theta) = \omega e^{-\omega \theta}$, $\theta \ge 0$, with average lifetime $\langle \theta \rangle = \omega^{-1}$. The mRNA age distribution is thus also exponential and can be written as

$$ \phi_a(a) = \omega e^{-\omega a}, \quad a \ge 0. \qquad (3) $$
In a biological cell, mRNA can be found in different states: newly synthesized mRNA empty of ribosomes, mRNA in a transient state evolving from an empty to stationary state, and finally mRNA in the stationary state of translation. We divide the whole population of mRNA into sub-populations, denoted k-somes containing exactly k ribosomes. The first four categories containing from one to four ribosomes, respectively, are called monosomes, disomes, trisomes and tetrasomes.

We will show in the next sections that these k-somes are more sensitive to mRNA degradation than the whole population of intracellular mRNAs (polysomes) and that these independent density profiles provide sufficient information for estimating all relative translation kinetic rates in a single experiment.

## Theoretical results

In this section we present an analysis of the k-somes within the framework of the ballistic model.

### Distribution of the k-some population

The number of ribosomes on a particular mRNA depends on its age a: if $a < T(L)$, ribosomes present on the mRNA at the time of the experiment have been loaded during the time interval a preceding the experiment while for $a > T(L)$, they have been loaded during the time interval T(L) before the experiment. The loading process being of the Poisson type, the probability that an mRNA of age a be loaded with k ribosomes is therefore given by

$$P(k|a) = \begin{cases} \frac{(\alpha a)^k}{k!} e^{-\alpha a}, & a \le T(L), \\ \frac{(\alpha T(L))^k}{k!} e^{-\alpha T(L)}, & a \ge T(L). \end{cases} \quad (4)$$

Using (3), the probability $P_k$ that an mRNA be a k-some irrespective of its age is then given by (see S1 Text for an explicit analytical expression)

$$P_k = \int_{0}^{\infty} P(k|a) \phi_a(a) da. \quad (5)$$

### k-some density profile along mRNA

Another important result about Poisson processes is that, knowing that k events have occurred during a given time interval, the probability density of their occurrence times is uniform on that time interval (see for example Ref. [16] for a rigorous statement of that theorem). Given the one-to-one correspondence between the ribosome initiation time and its location at the time of the Ribo-seq experiment, using the ribosome equation of motion, $dT(x)/dx = 1/p(x)$, the density profile $\rho_k(x|a)$ of a k-some knowing that its age is a can be written as

$$\rho_k(x|a) = \begin{cases} \frac{k}{ap(x)} H(x(a) - x), & a \le T(L), \\ \frac{k}{T(L)p(x)} , & a \ge T(L), \end{cases} \quad (6)$$

where $H(x)$ stands for the Heaviside function. In the first line of (6), $x(a)$ is the maximal distance that a ribosome can possibly travel along an mRNA with an age $a < T(L)$. It is given implicitly by $a = \int_{0}^{x(a)} dy/p(y)$. (6) is easily interpreted: for recently synthesized mRNA such that $a < T(L)$, a ribosome reaches at most the distance $x(a)$ along the mRNA. The density is
therefore zero beyond that limit (hence the Heaviside function) and the system behaves as if the mRNA had an effective length $x(a) < L$. Moreover, the ribosome density is proportional to the time ribosomes spend at location $x$ which is inversely proportional to their local velocity $p(x)$.

According to (6), it is possible to distinguish between two states for mRNA with translating ribosomes: a transient state, for which young mRNA with $a < T(L)$ have densities $\rho_k(x|a)$ that depend explicitly on the age $a$, and a stationary state for which $a > T(L)$ and $\rho_k(x|a)$ no longer depends explicitly on $a$. We will see later that the proportion of mRNA in one state versus the other is critical in characterizing the age-averaged density. Integrating any of these densities over the entire mRNA length yields, as expected, the number of ribosomes, $k$.

The density profile of a k-some irrespective of mRNA age—which is experimentally accessible for small values of $k$—is given by

$$ \rho_k(x) = \int_0^\infty \rho_k(x|a)p(a|k) da, \qquad (7) $$

where $p(a|k)$ is the age distribution among k-somes given by Bayes theorem, $p(a|k) = \phi_a(a)P(k|a)/P_k$. The analytical expression for $\rho_k(x)$ is somewhat complex and shall be provided and analyzed in section.

## Polysome density profile along mRNA

The ribosome density profile obtained for the entire mRNA population (polysomes) may be calculated as the average over $k$ of the k-some density profiles

$$ \rho(x) = \sum_{k=0}^{\infty} \rho_k(x) P_k = \frac{\alpha}{p(x)} e^{-\omega T(x)}. \qquad (8) $$

Integrating over $x$ leads to the average number of ribosomes $\langle k \rangle$ loaded on an mRNA

$$ \langle k \rangle = \frac{\alpha}{\omega} \left( 1 - e^{-\omega T(L)} \right). \qquad (9) $$

The two former equations, (8) and (9), are in agreement with Ref. [4], showing the consistency of our approach. Note that the result (9) can also be obtained from $\langle k \rangle = \sum_{k=0}^{\infty} kP_k$. Note also that, in the limit of infinite lifetime ($\omega = 0$), the polysome density (8) becomes $\rho^{\infty}(x) = \alpha/p(x)$, a result that simply expresses the conservation of the local particle current $j(x) = \rho(x)v(x) = \text{cst}$ in the stationary regime, and that can indeed be derived in the TASEP low density phase [17].

## Analysis of finite lifetime effects

In this section we rationalize the parametric dependence of the previous expressions in the relevant physical regimes. We shall use $T(L)$ as a time scale throughout the rest of the paper. This allows us to define dimensionless rates (with tildes) as original ones multiplied by $T(L)$

$$ \tilde{\alpha} = \alpha T(L); \quad \tilde{\omega} = \omega T(L); \quad \tilde{p}(x) = p(x)T(L). \qquad (10) $$

## Values of kinetic rates in various organisms

Since we are studying the effects of finite mRNA lifetime, we report in Table 1 some literature values for the kinetic parameters and in Table 2 their dimensionless counterparts that prove essential in determining typical biological regimes. We note that typical mRNA lifetimes
<table><caption>**Table 1. Translation kinetic rates and mRNA properties in bacteria, yeast and mammals.** Medians of gene length *L* are calculated from the entire genomes set available on the NCBI database. mRNA half-lives *t*<sub>1/2</sub> [18–21], elongation rates *p* [22–25] and initiation rates *α* [9, 26] are taken from the literature. *T*(*L*) is approximated by *L*/*p* (unit length is expressed in codons).</caption><thead><tr><th>Species</th><th>mRNA t<sub>1/2</sub> (median)</th><th>p (s<sup>-1</sup>)</th><th>L(codons)</th><th>α (s<sup>-1</sup>)</th><th>T(L) (s)</th></tr></thead><tbody><tr><td>E. coli</td><td>5 min</td><td>10–20</td><td>256</td><td>-</td><td>13–26</td></tr><tr><td>S. cerevisiae</td><td>20 min</td><td>3–10</td><td>399</td><td>0.1</td><td>40–133</td></tr><tr><td>M. musculus</td><td>7 h</td><td>6</td><td>493</td><td>-</td><td>82</td></tr><tr><td>H. sapiens</td><td>9 h</td><td>3</td><td>481</td><td>0.06</td><td>160</td></tr></tbody></table>

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

$t_{1/2} = 1/\omega$ for different species can span two orders of magnitude, while molecular quantities vary significantly less (e.g. a factor 7 for *p*, a factor 2 for *L*, a factor 10 for $T(L)$).

## Polysomes are not sensitive to mRNA degradation

In this subsection, we discuss the effect of degradation on the whole population of mRNA. In doing so, we shall define $f(\tilde{\omega})$ as the average number of ribosomes $\langle k \rangle$ normalized by the same quantity in the infinite lifetime limit, $\langle k \rangle_{\infty} = \tilde{\alpha}$ (see Eq (9)):

$$f(\tilde{\omega}) = \frac{\langle k \rangle}{\langle k \rangle_{\infty}} = \frac{\langle k \rangle}{\tilde{\alpha}} = \frac{1}{\tilde{\omega}}(1 - e^{-\tilde{\omega}}). \quad (11)$$

This function can be interpreted in terms of a birth and death process (see S2 Text). It only depends on $\tilde{\omega}$ and is plotted in Fig 2 in semilog scale.

As expected, $f(\tilde{\omega})$ decreases with increasing $\tilde{\omega}$: the average number of ribosomes per mRNA decreases as the mRNA lifetime becomes shorter. Nonetheless, biological values from the literature presented in Tables 1 and 2 show that in yeast and bacteria, $\tilde{\omega} \sim 10^{-2}$, while in mammals $\tilde{\omega} \sim 10^{-3}$. Therefore, $\tilde{\omega} \ll 1$ for the majority of genes and $f(\tilde{\omega}) \sim 1 - \tilde{\omega}/2$ remains very close to unity, showing that degradation has very little impact on $\langle k \rangle$, the most relevant parameter characterizing polysomes.

To strengthen even more this argument, let us discuss the polysome density profile, $\rho(x)$, which is one of the experimentally accessible quantities measured in Ribo-seq experiments, and compare it to its infinite lifetime limit, $\rho^{\infty}(x) = \alpha/p(x)$. We see from Eq (8) that the influence of mRNA degradation on $\rho(x)$ is maximal at the last site (i.e., $x = L$), which leads to the following bound

$$\forall x \in [0, L], \quad \frac{|\rho(x) - \rho^{\infty}(x)|}{\rho^{\infty}(x)} \le 1 - e^{-\tilde{\omega}} \approx \tilde{\omega}, \quad (12)$$

where the approximation holds for $\tilde{\omega} \ll 1$. For $\tilde{\omega} = 10^{-1}$, an upper bound of the values indicated in Table 2, replacing the polysome density by its infinite lifetime limit would lead to a

<table><caption>**Table 2. Calculated dimensionless parameters in bacteria, yeast and mammals from parameters taken from Table 1 (note that $\omega = 1/t_{1/2}$).** The quantity $R_1(0)$ will be defined in the next section, see Eq (21).</caption><thead><tr><th>Species</th><th>ω = αT(L)</th><th>α = αT(L)</th><th>ω/α</th><th>R<sub>1</sub>(0)</th></tr></thead><tbody><tr><td>E. coli</td><td>0.04–0.07</td><td>-</td><td>-</td><td>-</td></tr><tr><td>S. cerevisiae</td><td>0.02–0.08</td><td>4–13</td><td>6×10<sup>-3</sup></td><td>0.3–3700</td></tr><tr><td>M. musculus</td><td>0.002</td><td>-</td><td>-</td><td>-</td></tr><tr><td>H. sapiens</td><td>0.003</td><td>10</td><td>3×10<sup>-4</sup></td><td>5</td></tr></tbody></table>

https://doi.org/10.1371/journal.pcbi.1011522.t002
<figure><img src="image_5.png" /><figcaption>Fig 2. Horizontal semilog plot of the function $f(\tilde{\omega}) = \langle k \rangle / \tilde{\alpha}$ defined in Eq (11) vs $\tilde{\omega} = \omega T(L)$ (solid red line) and asymptotic behavior when $\tilde{\omega} \to 0$, $f(\tilde{\omega}) \sim 1 - e^{-\frac{\alpha}{2}}$ (dashed blue line). The function has been plotted over 5 biologically relevant decades: $\tilde{\omega} \in [10^{-5}, 1]$.</figcaption></figure>

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

relative error of no more than 10%. Therefore, no modifications due to mRNA degradation should be detectable in polysome Ribo-seq experiments. This analysis shows that previous works [2, 9] estimating the hopping rates in S. cerevisiae from polysomes with a model that does not take into account the mRNA finite lifetime provide a valid approximation. At variance with polysomes, we shall demonstrate in the next sections that, under the same conditions, k-somes are more sensitive to degradation than polysomes and constitute a key tool to determining kinetic rates.

### k-some analysis

**k-some distribution.** In terms of the dimensionless parameters $\tilde{\alpha}$ and $\tilde{\omega}$, the probability $P_k$ (see S1 Text) that an mRNA be a k-some is given by

$$P_k = \frac{\tilde{\omega}}{\tilde{\alpha} + \tilde{\omega}} \left( \frac{\tilde{\alpha}}{\tilde{\alpha} + \tilde{\omega}} \right)^k \frac{\gamma(k+1, \tilde{\alpha} + \tilde{\omega})}{k!} + \frac{\tilde{\alpha}^k}{k!} e^{-(\tilde{\alpha} + \tilde{\omega})}, \quad (13)$$

where details about the incomplete Gamma function, $\gamma(k, z)$, may be found in S3 Text. In Fig 3, this distribution is plotted for typical values of the parameters: $T(L) = 200s$, $L = 100$ codons, $\alpha = 0.06s^{-1}$ ($\tilde{\alpha} = 12$) and different values of $\tilde{\omega}$ between 0 (no degradation) and 1 (high degradation). A log-lin scale is used in the graph (A) to better observe the evolution of the k-some population with small $k$. As $\tilde{\omega}$ increases, we see that the populations of empty mRNA (zero-somes) and of the mono- through tetrasome increase quasi linearly with $\tilde{\omega}$ and that the
<figure><img src="image_6.png" /><figcaption>Fig 3. The k-some population distribution $P_k$ for finite values of $\tilde{\omega}$ (from 0 to 1) and other parameters fixed: $T(L) = 200s$, $L = 100$ codons and $\alpha = 0.06s^{-1}$. The value $\langle k \rangle = \tilde{\alpha} = 12$ corresponds to the maximum of the distribution at $\tilde{\omega} = 0$. (A) In the lin-log scale graph we can appreciate the increase of the population of empty mRNA ($k=0$) and small value of $k$ as $\tilde{\omega}$ increases. (B) Using the lin-lin scale, we can see that the two maxima (at $k = \tilde{\alpha}$ and $k = 0$) equate for $\tilde{\omega} = 0.8$. Expression (14) for $P_k$ provides an excellent approximation for $k \le 4$, except for the highest two values of $\tilde{\omega}$ (for which the approximation remains good with less than 20% error for $k \le 1$).</figcaption></figure>

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

distribution progressively displays a plateau for low values of $k$. This result is best understood by carrying out a linear expansion of $P_k$ in $\tilde{\omega}$ for $\tilde{\alpha} \gg k$:

$$P_k \approx P_k^\infty + \frac{\tilde{\omega}}{\tilde{\alpha}}, \quad \text{with} \quad P_k^\infty = \frac{\tilde{\alpha}^k e^{-\tilde{\alpha}}}{k!}. \qquad (14)$$

For values of $\tilde{\alpha} \sim 10$, the first term of this expansion, $P_k^\infty$—that is of the Poisson type, as expected given the nature of the ribosome loading process—remains very small compared to the linear term in $\tilde{\omega}$. We therefore have $P_k \sim \tilde{\omega}/\tilde{\alpha}$, a result both linear in $\tilde{\omega}$ and independent of $k$ that explains the observed plateau. For $\tilde{\omega} < 10^{-1}$, when $k$ increases, however, the k-some distribution becomes poissonian again and $P_k \sim P_k^\infty$, as we may observe in Fig 3.

In the graph (B), the same data are plotted on a lin-lin scale and with a broader range for $k$. For $\tilde{\omega} < 10^{-1}$, we observe a nearly symmetric Gaussian distribution that mainly corresponds to the Poisson term $P_k^\infty$. It is centered around a maximum close to $k = \langle k \rangle_\infty = \tilde{\alpha} = 12$, the value of the average number of ribosomes in the infinite lifetime limit. As expected, most mRNAs are then loaded with that number of ribosomes when $\tilde{\omega} < 10^{-1}$. For higher values of $\tilde{\omega}$, however, typically when $\tilde{\omega} > \tilde{\omega}_{inv}$ where $\tilde{\omega}_{inv} = \ln(\sqrt{\tilde{\alpha}/2\pi} + \frac{1}{2})$ (see S1 Text), an inversion of the populations occurs: the proportion of small order k-somes takes over the population of those for which $k \sim \tilde{\alpha}$ (local maximum). This result corresponds to a regime where mRNAs do not live long enough for ribosomes to reach the exit, hence reducing their number and favoring k-some populations with small $k$. For biologically relevant values ($\tilde{\alpha} \sim 10$), $\tilde{\omega}_{inv}$ is on
the order of unity, in good agreement with the value $\tilde{\omega} \sim 1$ for which the polysome density is significantly impacted by degradation (see (12)).

**k-some densities.** One might suspect from the previous results that, due to their statistical prominence when $\tilde{\omega} \sim 1$, mono- through tetrasomes (k-somes with small k) might play a decisive role in the determination of the kinetic parameters only once that condition is fulfilled. This is, however, not true for, as we shall see, k-some densities become sensitive to mRNA degradation for smaller values of $\tilde{\omega}$ than those ensuring their dominance. Therefore, inasmuch as a statistically relevant sample of them (say thousands) is available in Ribo-seq data, densities may be analysed to determine the kinetic parameters. Hundreds of millions of mRNA sequences are typically involved in Ribo-seq experiments. This means that for $\tilde{\alpha} \sim 10$ and as soon as $\tilde{\omega} > 10^{-4}$, thousands of monosomes and even more of k-somes with $k > 1$ are available for analysis, enough to reach a good statistical accuracy.

In terms of dimensionless parameters, the exact k-some density (7) reads

$$ \rho_k(x) = \frac{\tilde{\alpha}^k e^{-(\tilde{\alpha}+\tilde{\omega})}}{P_k \tilde{p}(x) (k-1)!} + \frac{\tilde{\omega} \left(\frac{\tilde{\alpha}}{\tilde{\alpha}+\tilde{\omega}}\right)^k \gamma(k, \tilde{\alpha}+\tilde{\omega}) - \gamma(k, (\tilde{\alpha}+\tilde{\omega})\tau(x))}{P_k \tilde{p}(x) (k-1)!} \quad (15) $$

where $\tau(x) = T(x)/T(L) \in [0, 1]$. This expression is composed of two terms, each corresponding to those obtained in (6). The first one represents the contribution to the k-some density of mRNAs that have reached their stationary state of translation in the sense that a ribosome may have had the chance to exit them prior to the experiment. We shall refer to that term as $S_k(x)$. The second term corresponds to mRNAs that are still in their transient state of translation: the time elapsed between their synthesis and the Ribo-seq experiment is too short for a ribosome to reach their exit codon. We shall call that term $F_k(x)$.

To quantify the influence of mRNA degradation on k-somes, we shall define a transient to stationary mRNA density ratio as $R_k(x) = F_k(x)/S_k(x)$. A mathematical analysis of this ratio is done in S4 Text. Its maximal value, achieved at $x = 0$, is given by

$$ R_k(0) = \tilde{\omega} \frac{\gamma(k, \tilde{\alpha} + \tilde{\alpha})}{(\tilde{\omega} + \tilde{\alpha})^k} e^{\tilde{\omega} + \tilde{\alpha}} \approx \frac{\tilde{\omega}}{k P_k^\infty}, \quad (16) $$

where $P_k^\infty$ is given in (14). The last approximation is valid for small values of $\tilde{\omega}$ and $\tilde{\alpha} \gg k$ but the exact expression, valid for all values of $\tilde{\alpha}$ and $\tilde{\omega}$, allows us to analyse the impact of degradation in the entire parameter space. The fact that $R_k(x)$ achieves a maximal value at $x = 0$ means that k-some densities are more sensitive to degradation near the entrance site than towards the end, at variance with polysomes for which it is the opposite (see Section). We shall see below that the approximate value of $R_k(0)$ appears as a natural parameter in the comparison of the k-some density to its infinite lifetime limit.

To illustrate the global trends of k-some densities without further complications, and as a non uniform elongation rate would just modulate the k-some density without affecting its global behavior, we shall now assume a constant (codon independent) elongation rate, $\tilde{p}(x) = \tilde{p}$. In that case, the dimensionless time is simply $\tau(x) = x/L$ and the sole term of (15) that depends on $x$ is $F_k(x)$. This term is responsible for the global exponential decay of the densities observed in Fig 4 where $\rho_k(x)$ has been plotted for $k = 1$ to 4. In that figure, density profiles in solid lines are given by (15) (the discrete version of which can be found in S5 Text) with $L = 100$ codons, $\alpha = 0.06s^{-1}$, $p = 0.5s^{-1}$ and thus $T(L) = 200s$. Dimensionless parameters are thus $\tilde{\alpha} = 12$ and $\tilde{p} = 100$.

In each graph, dashed lines represent the density profiles for an infinite mRNA lifetime ($\tilde{\omega} = 0$) which are flat with magnitude $\rho_k^\infty = k/L$. In that same limit, the polysome density is also flat, with magnitude $\rho^\infty = \alpha/p$. As an important consequence, when $\tilde{\omega} = 0$, all k-some
<figure><img src="image_7.png" /><figcaption>Fig 4. Plot of the *k*-some densities ρ_k for *k* = 1 to 4 versus the genomic coordinate *x* obtained from (15): Monosome (blue), disome (orange), trisome (green) and tetrasome (red). Black curves correspond to polysome densities. Parameters are α = 0.06s⁻¹, *p* = 1/2s⁻¹ and *L* = 100 codons (π = 12). Solid lines are obtained from (15) for increasing finite values of (A) ω = 5.10⁻⁷s⁻¹ (ω̃ = 10⁻⁴), (B) ω = 5.10⁻⁵s⁻¹ (ω̃ = 10⁻²) and (C) ω = 5.10⁻³s⁻¹ (ω̃ = 1); while dashed lines in each panel are obtained for ω = 0s⁻¹ (infinite mRNA lifetime).</figcaption></figure>

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

and the polysome profiles are linearly related. We shall see in section that this precludes the
determination of all kinetic parameters in that limit.

From panel (A) to (C), the values of tilde{ω} used in the successive panels are spaced by a factor
of two orders of magnitude: (A) ω̃ = 10⁻⁴, (B) ω̃ = 10⁻² and (C) ω̃ = 1 corresponding to
555h, 5h and 3 min lifetimes, respectively. Below, we briefly describe the main features of the
densities ρ_k(x). A more detailed analysis is provided in S6 Text.

In Fig 4(A), all density profiles are close to the ̃ω = 0 limit and we barely notice a small
deviation for monosomes, at the start of the filament. Thus, in this regime of low degradation,
we do not expect *k*-somes to bring more information than the polysomes in a usual Ribo-seq
experiment.

In Fig 4(B), we do not see any effect on the polysome profile, but do see a clear effect on the
*k*-somes. Note that, as *k*-some profiles are no longer linearly related, a fit of the kinetic rates is
now possible. This decoupling indeed provides extra information for determining the
unknown model parameters. Importantly, this intermediate regime (̃ω ~ 10⁻²) falls into the
range of biological values. Additionally, we observe that the ordering of *k*-some densities
change with the position along the mRNA: while densities obey $ρ_1 > ρ_2 > ... > ρ_4$ close to the initiation site, they switch in reverse order after a few sites.

In Fig 4(C), degradation is high ($\tilde{\omega} = 1$). All density profiles are modified, even the polysomes. This regime is probably not relevant for typical translation processes, but we describe it for the sake completeness and because it might be relevant under extreme cellular conditions.

Let us now turn to the sensitivity of k-some density profiles to mRNA degradation. Taking advantage of the smallness of $\tilde{\omega}$ compared to other parameters, we carry out a linear expansion of $ρ_k(x)$ in $\tilde{\omega}$. Yet, as we may observe in Fig 4, k-some densities are more sensitive to degradation near the ribosome initiation site ($x = 0$) than further along the mRNA. More details about the subtle question of the position dependent sensitivity of k-some densities are provided in S6 Text, where the intersection of the k-some density profiles with their infinite lifetime limit (flat profiles) is analysed. Hereafter, we evaluate their sensitivity to $\tilde{\omega}$ by calculating their relative difference with their infinite lifetime limit at $x = 0$. For $\tilde{\alpha} \gg k$, we obtain

$$ \frac{|\rho_k(0) - \rho_k^\infty(0)|}{\rho_k^\infty(0)} \approx \frac{\tilde{\omega}}{kP_k^\infty} \approx R_k(0). \quad (17) $$

Thus, in the biologically relevant limit of small $\tilde{\omega}$ values, the transient to stationary mRNA ratio $R_k(0)$ controls the impact of degradation on k-some densities near the initiation site.

In conclusion, the ballistic model mainly displays three regimes of degradation: (i) low degradation (LD), where k-somes and polysomes are linearly related and where k-somes are therefore not expected to provide more information than polysomes; (ii) an intermediate regime (ID), into which some genes are likely to fall, where k-some densities are already sensitive to degradation and no longer linearly related while polysomes are still unaffected; (iii) a regime of high degradation (HD), not relevant for typical translation processes, in which all profiles are strongly modified by degradation, polysomes included.

## Comparison with experiments

### Ribo-seq experiments with k-somes

In this section, we compare the phenomenology of the ballistic model for k-somes to experimental data for a gene coding for histones. This comparison will show that the effect of degradation is visible in certain genes and the sensitivity to degradation is enhanced for k-somes with respect to polysomes, thus justifying the introduction of this new experimental setup.

Fig 5 displays in panel (A) the polysome fractioning procedure and in panels (B-C) the k-some densities, $ρ_k$ (for $k = 1$ to 4), versus the genomic coordinate obtained from Ribo-seq profiling experiments for the histone gene HIST1H2BM (symbols) and from a fit of these data using the homogeneous ballistic model (solid lines).

The polysomes were fractionated following the method of [27] in order to separate the mRNAs bound either to one, two, three or four ribosomes. This yields four subsets of mRNA respectively termed mono-, di-, tri-, and tetra-somes (Fig 5(A)). For each fraction, we extracted all mRNA segments covered by a ribosome (called Ribosome Protected Fragments (RPF)), and prepared a library for sequencing as described in [28]. Then short read sequencing was performed as for RNA-seq using an Illumina platform. The sequencing data for each library (mono-, di-, tri-, and tetra-somes), which contains raw sequencing reads, are processed the same way, one independently of the others (using a method adapted from [29]). For a given mRNA/gene, the goal of bioinformatic processing is to tabulate how many reads cover a given position in the mRNA coding sequence. After cleaning, each read sequence is aligned to the mRNA sequence to determine the position of origin of the RPF—a classical task of sequence comparison called mapping—using efficient software [30]. Only reads having a high
<figure><img src="image_8.png" /><figcaption>**Fig 5.** (A) Absorbance profile of a representative experiment. RNA populations were separated on a 15–50% sucrose gradient under 35,000rpm centrifugation for 3.5h. Free RNA, 40S ribosomal subunit, 60S ribosomal subunit, mRNA charged with 1,2,3...ribosomes were separated through a 254 nm UV spectrometer and collected with an ISCO fractionation system. Fractions corresponding to mRNA charged with 1,2,3 or 4 ribosomes were pooled and subjected to RPF extraction and cDNA construction. (B) *k*-some densities ρₖ (k = 1 to 4) obtained from the homogeneous ballistic model (solid lines) and from Ribo-seq profiling experiments (symbols) are plotted versus the genomic coordinate (in codon unit). The polysome density is plotted in black for the ballistic model only (lack of experimental data). Experimental data correspond to the histone gene HIST1H2BM and the fit is based on the ballistic model (see text for details) with fixed parameters, α = 0.06s⁻¹ and ω = 1/60 min⁻¹ ≈ 3.10⁻⁴ s⁻¹. The fit then yields T(L) = 178s and hence p = 0.7 codons/s. (C) Same *k*-some densities are normalized to 1, i.e. ρ̄ₖ = ρₖ/k, in order to compare the discrepancies with each-other and show that we retrieve the characteristic crossing between ballistic *k*-somes curves.</figcaption></figure>

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

quality alignment and a unique origin location are further processed. One records the match-
ing position of the first base of the read. On the first position of each read, we applied a shift of
12 nucleotides to assign the read to the P-site position of ribosomes [29]. It has been observed
from previous analyses that 12 nucleotides is by far the most frequent distance between the
read start and the P-site [31], which determines the position of the codon being translated.
Finally, we record the counts for each P-site covered position of the mRNA in a table.

In the panel (B) of Fig 5, Ribo-seq data (cross symbols) have been subsequently smoothed
over 19 codons by averaging the signal within a sliding bin and normalized such that the integral
of each k-some density ρk corresponds to k. As the genomic coordinate increases, we observe a
global exponential-like decay of the densities, a phenomenon predicted by the ballistic model
when mRNA degradation is taken into account (see Fig 4). Furthermore, another clear signature
of mRNA degradation may be observed in the panel (C) of Fig 5 where we display the Ribo-seq
k-some densities *normalized to unity*. These profiles show clearly that near the entrance of the
mRNA (over a genomic distance of 30–40 codons), k-some density profiles decrease with
increasing *k* while they increase with *k* beyond that distance. This crossing of the *k*-some density
curves is another prediction of the ballistic model in the presence of mRNA degradation.
Encouraged by the similarities between k-some Ribo-seq data and the corresponding ballistic trends noted above, we shall now try to compare them more quantitatively. We shall do so by restricting this first comparison to the homogeneous version of the ballistic model, i.e., $p(x)$ constant independent of the codon, and to smoothed experimental data. This choice is partly dictated by its simplicity of implementation, but also and foremost, by the lack of polysome data in the Ribo-seq set available for this histone gene (polysome data prove to be crucial in determining the codon-dependent elongation rates, $p_i$, as we shall explain in the next section). In their absence we choose to fit an average (and thus constant) elongation rate $p$ while fixing the two other parameters, $\alpha$ and $\omega$, to reasonable values found in the literature. To do so, we use the specific half-life of histone mRNA ($\tau_{1/2}$ = 40 min) to obtain (by definition) $\omega = \ln(2)/\tau_{1/2} \approx 1/60 \text{ min}^{-1} \approx 3.10^{-4} \text{ s}^{-1}$ [32] and use the initiation rate for *Homo sapiens* given in [26]: $\alpha = 0.06s^{-1}$. Finally, $p$ is estimated through $\mathcal{T}(L) = L/p$ in order to optimize the fit for all four studied values of $k$. We find $\mathcal{T}(L) = 178$ s, a value that leads to an average number of ribosomes on the mRNA in the stationary regime $\langle k \rangle = \alpha\mathcal{T}(L) = 0.06 \times 178 \approx 10$. This value is the average of $\mathcal{T}(L)$ values estimated for specific values of $k$: $\mathcal{T}(L) = 147$s for monosomes, 168$s$ for disomes, 190$s$ for trisomes and 206$s$ for tetrasomes. Given the crudeness of this first approach where the codon-dependency of the elongation rates is not taken into account, the dispersion of these results seems reasonable.

Overall, we observe that the qualitative trends exhibited by the Ribo-seq data are described by the ballistic model and that the inversion of the *k*-some density profiles near the start of the mRNA gives a clear signature that translation in these experiments is influenced by mRNA degradation. Finally, in contrast to the *k*-somes, the theoretical density for polysomes (black solid line in the intermediate panel of Fig 5), which has no experimental counterpart here, is nearly flat and barely influenced by degradation: this is in agreement with the regime of intermediate degradation presented in Fig 4(B). This result provides strong further motivation for the need to split the population of mRNA into *k*-somes in order to increase the sensitivity to mRNA degradation. Note that from the analysis of $f(\tilde{\omega})$ done in section, it was expected that in the present case mRNA degradation would have little influence on the polysomes density since $\tilde{\omega} = 3.10^{-4} \times 178 \approx 0.05 \ll 1$. Finally, we conclude this section by remarking the need to introduce an inhomogeneous $p(x)$ to explain the strong and systematic fluctuations along the profiles that cannot be accounted for by experimental fluctuations. In the absence of polysome data to complete the *k*-some set, we deem it premature to attempt such a fit here (experiments in progress).

# Fitting the kinetic rates from ribosome density profiles

## Limitations in using solely polysomes to fit kinetic rates

Usual Ribo-seq experiments involve polysomes only and their modeling does not account for the mRNA degradation rate, $\omega$. This is justified inasmuch as, for most living cells, the polysome density is barely sensitive to degradation and this approximation does not lead to significant errors in the estimation of the elongation rates. Yet, the main problem in estimating the kinetic rates from polysomes alone is that there are more unknowns than equations (or equivalently, more model parameters than experimentally measured quantities). Considering henceforth the inherent discreteness of the mRNA, we shall now use $p_i$ as the elongation rate at site (codon) $i \in \{1, \dots, L\}$ and $\rho_i$ for the polysome density at that site instead of their continuous version. Now, even if we assumed the degradation rate $\omega$ to be known from independent experiments, we would still need to estimate $(L+1)$ parameters: the set of $\{p_i\}$ and the initiation rate $\alpha$. As the polysome density $\rho_i$, on the other hand, only provides $L$ data, we are left with an extra unknown.
Another problem in estimating the kinetic parameters is that the number of mRNAs from which the Ribo-seq profile is obtained is not known, which prevents a proper normalization of the Ribo-seq density. A parametric estimation must therefore rely on a density $\bar{\rho}(x)$ normalized to unity, $\int_0^L \bar{\rho}(x)dx = 1$, instead of $\langle k \rangle$. This density is given by

$$ \bar{\rho}(x) = \frac{\omega e^{-\omega \mathcal{T}(x)}}{p(x) (1 - e^{-\omega \mathcal{T}(L)})} \quad (18) $$

Although this density no longer depends on $\alpha$, this does not solve the problem, since the new normalisation constraint, whose discrete version is $\sum_{i=1}^L \bar{\rho}(i) = 1$, also removes one equation for the input parameters. A strategy found in the literature to circumvent this problem is to normalize Ribo-seq profiles by a mean ribosome velocity measured from independent experiments [2]. This method suffers from several drawbacks, however: (i) the ribosome velocity that is used is averaged over an entire genome instead of being gene specific, (ii) results are typically obtained for species different from the one involved in the Ribo-seq experiment and (iii) experiments are done independently in different conditions. Altogether, the normalization obtained from this procedure is expected to be subject to large uncertainties.

## k-some improved fit of kinetic rates

Half-lives ($\tau_{1/2}$) of individual mRNAs have been successfully measured in mouse embryonic stem cells [20]. By applying the same method, any mRNA half-life can be measured as well. We shall henceforth consider $\omega = \ln(2)/\tau_{1/2}$ as a known kinetic time scale. The normalised polysome density does not depend on $\alpha$ (Eq 18), but on the $p_i$'s only and, as we have shown, it is not strongly affected by mRNA degradation. It therefore constitutes a suitable quantity to determine the variation along the mRNA of the elongation rates, $p_i$, independently of the other parameters. Furthermore, in section, we have shown that a k-some Ribo-seq data analysis provides additional information on the kinetic parameters, especially when the mRNA lifetime decreases for, in that case, k-some densities become less and less linearly related. Yet, when the mRNA degradation becomes too important (large $\omega$), k-some density profiles vanish near the end of the mRNA rendering their analysis quite difficult. The same is true of polysome densities. This, however, happens for parameters far outside the biological range.

In a first approach, we therefore decide to combine the information provided by both the polysome normalized density and the monosome density since among k-some profiles, it is the most sensitive to the mRNA lifetime (see section and S4 Text).

## Optimizing the fit with respect to the degradation rate

We now implement the strategy outlined above by detailing how to fit the initiation rate $\alpha$ and the elongation rates $p_i$ of all the codons constituting the mRNA. To find these kinetic parameters, we minimize the quadratic difference between the ballistic model predictions and, simultaneously, the Ribo-seq data for the normalized (to one) polysome density and the monosome density. To be more specific, we minimise the quantity

$$ \mathcal{X}^2 = \chi_{\text{polysome}}^2 + \chi_{\text{1-some}}^2, \quad (19) $$

where $\chi_{\text{polysome}}^2$ is the sum of the relative quadratic errors of the normalized polysome density on each codon:

$$ \chi_{\text{polysome}}^2 = \sum_{i=1}^{L} \left| \frac{\rho_i^{\text{data}} - \rho_i^{\text{ballistic}}}{\rho_i^{\text{data}}} \right|^2, \quad (20) $$
and where $\chi_{1-some}^2$ has a similar form for monosomes. The minimisation, performed by the L-BFGS algorithm (an adaptation of the Broyden-Fletcher-Goldfarb-Shanno quasi-Newton method for using a limited amount of computer memory) encoded in the Python 3 library Scipy [33], yields at once $\alpha$ and the L parameters $p_i$.

As we shall see, the accuracy of our fitting method improves when the values of the ratio

$$ \mathcal{R}_1(0) = \tilde{\omega}(e^{\tilde{\omega}+\tilde{\alpha}} - 1)/(\tilde{\omega} + \tilde{\alpha}), \qquad (21) $$

(see (16) and S4 Text) increase. As a first feasibility test, the minimisation procedure is performed on artificial data created by the ballistic model itself. The goal is to measure its accuracy in recovering the known input parameters as a function of $\mathcal{R}_1(0)$. The relative error between a parameter value found by minimisation, $Y^*$, and its exact value, $Y$, is defined by

$$ \text{Err}_Y = \frac{|Y - Y^*|}{Y}. \qquad (22) $$

We call it a “score function”. Beside $\text{Err}_\alpha$, and to avoid the wealth of score functions related to every single elongation rate $p_i$, we shall estimate two quantities: an average score function $\text{Err}_{\text{mean}} = \sum_i \text{Err}_{p_i}/L$, and the maximum of them all, $\text{Err}_{\text{max}} = \max_i \text{Err}_{p_i}$.

To test our minimisation procedure, we fix $\alpha = 0.08 \text{ s}^{-1}$. We then generate ten sets of L random $p_i$ values drawn from a uniform distribution in the range [0.2, 4] s$^{-1}$. Each set has its own $\mathcal{T}(L) = \sum_i p_i^{-1}$ value and $\omega$ is adjusted in such a way that the ratio $\mathcal{R}_1(0)$ has a desired value (here, an integer from 1 to 15). Now, to each set of $\alpha, p_i$ and $\omega$ values correspond ballistic polysome and monosome densities. These densities are the artificial data processed by minimisation of the $\mathbb{X}^2$ function (19) as explained above. This procedure yields the parameter values $\alpha^*$ and $p_i^*$ from which the score functions $\text{Err}_\alpha, \text{Err}_{\text{mean}}$ and $\text{Err}_{\text{max}}$ are then evaluated to check the accuracy of the method.

Fig 6 displays the 3 score functions for the ten density profiles sharing the same $\alpha = 0.08 \text{ s}^{-1}$ and the same $\mathcal{R}_1(0)$ (pale colors). The observed variability stems from the randomness of the $p_i$ sets. To smooth out this variability, score functions have subsequently been averaged over the ten random profiles (vivid colors). We can see that relative errors decrease when $\mathcal{R}_1(0)$ increases: they are less than 10% for $\mathcal{R}_1(0) \ge 3$ and remain nearly constant (around 1%) for $\mathcal{R}_1(0) \ge 10$. This confirms the relevance of $\mathcal{R}_1(0)$ as a good criterion for the fit quality, albeit in a suitable parameter range. In the presented fits, $\mathcal{R}_1(0) = 3$ corresponds to the parameters $\alpha = 0.08 \text{ s}^{-1}$, $\mathcal{T}(L) \approx 81 \text{ s}$ and $\omega \approx 5.10^{-4} \text{ s}^{-1} \approx \frac{1}{33} \text{ min}^{-1}$, which are all biologically reasonable rates for mRNA translation. This analysis therefore allows us to assess the feasibility of our proposed fitting method to determine the kinetic parameters of mRNA translation for a particular gene. We plan, in a forthcoming paper, to compare quantitatively the fitted values obtained with other previous methods using solely the polysome profile. To carry out this comparison thoroughly, we have planned to experimentally measure both the polysome and the k-some profiles with the same cell line, and compare the values obtained with these two approaches.

### Translation phase diagram: The effect of finite mRNA lifetime

We summarize our main results in Fig 7. This figure can be considered as a *phase diagram* of the ballistic model in the $(\tilde{\alpha}, \tilde{\omega})$, or (initiation, degradation), plane. It indicates the three phases already defined in section: low degradation (LD), where neither polysome nor k-some ribosome density profiles are affected by the finite mRNA lifetime, intermediate degradation (ID), where only k-some densities (with small $k$) are affected, and high degradation (HD), where both polysome and k-some densities are impacted. The LD and ID regimes are separated by the curve $\mathcal{R}_1(0) = 5$, above which k-somes feel a pronounced effect of degradation.
<figure><img src="image_9.png" /><figcaption>Fig 6. Score functions (22) of the fit with the ballistic model when $\alpha = 0.08s^{-1}$ and for 10 ballistic density profiles obtained by a random drawing from a uniform distribution of the elongation rates $p_i \in [0.2, 4]$ s$^{-1}$ (pale colors). Results averaged over the 10 profiles are displayed in vivid colors. $R_1(0)$ (21), is increased by increasing $\omega$ and $\mathcal{T}(L) \approx 81.3$ s on average.</figcaption></figure>

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

<figure><img src="image_10.png" /><figcaption>Fig 7. Phase diagram of the ballistic model in the plane $(\tilde{\alpha}, \tilde{\omega})$. The phase diagram is split in three degradation regimes with the lines of equation $\tilde{\omega} = 1 \sim \mathcal{O}(1)$ and $R_1(0) = 5 \sim \mathcal{O}(1)$. Biological values obtained for the histone HIST1H2BM (human) are given by the homogeneous fitting from monosome to tetrasome (Fig 5).</figcaption></figure>

https://doi.org/10.1371/journal.pcbi.1011522.g007
The ID and HD regimes are separated by the line $\tilde{\omega} = 1$, where polysomes are strongly affected by degradation, which becomes visible on the constant $\langle k \rangle$ curves in the ID phase and therefore on polysomes when $\tilde{\omega} > 0.1$, see also Fig 2.

The ID phase also turns out to be the most efficient regime for our fitting method. In this regime, the equations for the k-some densities are decoupled and the impact of degradation is not yet strong enough for k-some densities to be too low at the end of the mRNA transcript to obtain information about $p_i$ variations there. We have noted, however, that the quality of the fits doesn't decrease until $\tilde{\omega} = 5$ is reached, where the polysome density becomes of the order of $10^{-5}$ near the end of the mRNA and information on the corresponding elongation rates become undetectable.

The dashed lines represent curves of constant $\langle k \rangle$ (Eq (9)) in the $(\tilde{\alpha}, \tilde{\omega})$ plane. As degradation increases, we observe a departure from the value $\langle k \rangle_{\infty} = \tilde{\alpha}$ that is a good approximation in the intermediate regime. We also see that this value does not evolve substantially until the boundary between the ID and HD phases is approached. We observe that the line $\langle k \rangle = 1$ intersects the curves $R_1(0) = 5$ and $\tilde{\omega} = 1$ at the same point. An average value $\langle k \rangle = 1$ is then at the limit of feasibility for our proposed fitting method.

We have evidence from the results shown in Fig 7 that at least some histonic genes (green rectangle), as well as the average yeast genome value (S. cerevisæ, yellow rectangle), fall in the ID phase. The average for the *H. sapiens* genome, however, appears to fall close to the boundary between the LD and ID phase (red line) and may be at the limit of feasibility of our fitting method. Our analysis does suggest, however, other routes to determining the model parameters, e.g. by fitting the k-some repartition $P_k$, see Eq (13), an experimentally accessible quantity.

We postpone to a forthcoming paper a more detailed comparison with specific sets of genes. We are currently evaluating mRNA degradation rates by qPCR after 24h growth of the same synchronized SW620 cells with nocodazole followed by a time course experiment in the presence of actinomycin D to block transcription. Subsequently, we will perform RNA Trizol extraction, DNaseI treatment, M-MLV reverse transcription and qPCR. The decay rate will be obtained from the slope of a semilogarithmic plot of mRNA level (qPCR values) as a function of time.

It is interesting to note that in the ID regime, the method we propose to extract the kinetic parameters of translation exploits a monosome density profile that is experimentally accessible despite the possibly much lower probability of finding a monosome than a k-some with $k \approx \langle k \rangle_{\infty} \approx \tilde{\alpha} (P_1 \ll P_{\tilde{\alpha}}^{\infty})$. In that sense, a “deep” sequencing (billions of mRNAs corresponding to millions of cells) is necessary to study monosomes as thousands of them are needed in order to ensure a statistically relevant analysis (see discussion in section).

## Discussion and conclusion

In this paper, we present an experimental setup where mRNAs are separated, before being sequenced, according to the number of ribosomes attached to them. In this context, and for the first time to our knowledge, we develop a ballistic model of ribosome transport tailored to the analysis of such Ribo-seq experiments in the presence of mRNA degradation.

By investigating the mean ribosome number $\langle k \rangle$, as well as polysome and k-some densities, we characterize the system by identifying three regimes that emerge from a parametric analysis of the ballistic model. Provided that $\langle k \rangle_{\infty} = \tilde{\alpha} > 1$, we arrive at the following global picture: in the low degradation (LD) regime the effect of the finite mRNA lifetime on both $\langle k \rangle$ and the densities is negligible; in the intermediate degradation (ID) regime, the finite mRNA lifetime has a strong impact on k-some densities, but not on polysome densities (and therefore little impact on $\langle k \rangle$); in the high degradation (HD) regime, all densities, as well as $\langle k \rangle$, are affected
by the finite mRNA lifetime. The HD regime is delimited by $\tilde{\omega} = 1$ for which the time for a ribosome to cross the mRNA, $\mathcal{T}(L)$, is equal to the inverse of the degradation rate $\omega$. This limit probably does not correspond to any usual biological regime, which leads us to the conclusion that polysome densities should probably never be affected by finite mRNA lifetime (except perhaps in rare extreme conditions).

Given that the ID regime is the biologically relevant one, we can base a viable fitting strategy on the monosome and normalized polysome densities alone. Because the finite lifetime strongly affects the monosome density, and only slightly the polysome one, the two densities give access to complementary information. Whereas the monosome density exhibits the typical exponentially decreasing profile, depending on $\alpha$, $\omega$ and $\mathcal{T}(L)$, expected from a finite mRNA lifetime, the normalized polysome density is mainly sensitive to $p_i$ variations. More precisely, when normalized, the polysome density no longer depends on $\alpha$, and since $\omega$ is measured by another experiment, the only unknowns are the $p_i$. Furthermore, we have identified a quantity directly related to the effectiveness of the method, namely the ratio of transient to stationary monosome densities, $R_1(0)$, (21). Since the experiments presented in this article involve measuring the density profiles for both the polysomes and the k-somes with $k = 1, \dots, 4$, the model predictions for the disomes, trisomes, and tetrasomes (not used in the fitting procedure) provide a test of the consistency of the modeling approach we propose in the ID regime.

Using our ballistic model, we have analyzed preliminary Ribo-seq data for a human histonic gene. It is clear from this analysis that the Ribo-seq data display the same phenomenology predicted by the ballistic model in the ID regime: a decrease of the k-some Ribo-seq profiles at the end of the mRNA due to the degradation of mRNA in the late stages of translation. Thus we show for the first time the effect of degradation on Ribo-seq data.

In Fig 7, we present a phase diagram for the ballistic model in the $(\tilde{\alpha}, \tilde{\omega})$ plane as a global summary of the main results of our finite lifetime effect study and fitting method. This phase diagram displays the three regimes predicted by our analysis—LD, ID and HD—with the ID regime being considered the biologically relevant one, as well as the regime where our fitting method is effective.

In the modeling of real experimental data it will also be important to evaluate the influence of experimental noise on the reliability of our fitting method. It would be interesting to look at the effect of degradation with finite resources of ribosomes [34] and finite diffusion of ribosomes [35, 36]. Our method could also be applied to other phenomena where degradation potentially plays an important role, such as its coupling to codon bias [37], mRNA coverage [38], or ribosome drop off [39].

## Supporting information

S1 Fig. Quantifying the effect of excluded volume interaction in the ballistic model.
(PDF)

S1 Text. Analytical expression and analysis for $P_k$.
(PDF)

S2 Text. Function $f(\tilde{\omega})$: A birth and death approach.
(PDF)

S3 Text. Gamma function: Exact formulas and approximations.
(PDF)

S4 Text. Transient to stationary ratio for k-some densities.
(PDF)
**S5 Text. Discrete approach of the ballistic model.**

(PDF)

**S6 Text. Detailed analysis of the *k*-some densities.**

(PDF)

## Author Contributions

**Conceptualization:** Carole Chevalier, Jérôme Dorignac, Yahaya Ibrahim, Frédéric Geniet, Nils-Ole Walliser, John Palmeri, Andrea Parmeggiani, Jean-Charles Walter.

**Formal analysis:** Julie Ripoll, Eric Rivals.

**Investigation:** Armelle Choquet, Alexandre David.

**Supervision:** Andrea Parmeggiani, Jean-Charles Walter.

**Writing – original draft:** Carole Chevalier, Jérôme Dorignac.

**Writing – review & editing:** Jean-Charles Walter.

## References

1. Alberts BM, Johnson A, Lewis J, Morgan DO, Raff MC, Roberts K, et al. Molecular biology of the cell. 6th ed. New York [etc.]: Garland Science, Taylor & Francis Group; 2015.

2. Duc KD, Song YS. The impact of ribosomal interference, codon usage, and exit tunnel interactions on translation elongation rate variation. PLOS Genetics. 2018; 14(1):e1007166. https://doi.org/10.1371/journal.pgen.1007620 PMID: 30142215

3. Brar GA, Weissman JS. Ribosome profiling reveals the what, when, where and how of protein synthesis. Nature reviews Molecular cell biology. 2015; 16(11):651–664. https://doi.org/10.1038/nrm4069 PMID: 26465719

4. Valleriani A, Ignatova Z, Nagar A, Lipowsky R. Turnover of messenger RNA: Polysome statistics beyond the steady state. EPL (Europhysics Letters). 2010; 89(5):58003. https://doi.org/10.1209/0295-5075/89/58003

5. Siwiak M, Zielenkiewicz P. A Comprehensive, Quantitative, and Genome-Wide Model of Translation. PLOS Computational Biology. 2010; 6(7):e1000865. https://doi.org/10.1371/journal.pcbi.1000865 PMID: 20686685

6. Deneke C, Lipowsky R, Valleriani A. Effect of ribosome shielding on mRNA stability. Physical Biology. 2013; 10(4):046008. https://doi.org/10.1088/1478-3975/10/4/046008 PMID: 23883670

7. Reuveni S, Meilijson I, Kupiec M, Ruppin E, Tuller T. Genome-scale analysis of translation elongation with a ribosome flow model. PLoS computational biology. 2011; 7(9):e1002127. https://doi.org/10.1371/journal.pcbi.1002127 PMID: 21909250

8. Fernandes LD, Moura APSd, Ciandrini L. Gene length as a regulator for ribosome recruitment and protein synthesis: theoretical insights. Scientific Reports. 2017; 7(1):17409. https://doi.org/10.1038/s41598-017-17618-1 PMID: 29234048

9. Ciandrini L, Stansfield I, Romano MC. Ribosome Traffic on mRNAs Maps to Gene Ontology: Genome-wide Quantification of Translation Initiation Rates and Polysome Size Regulation. PLOS Computational Biology. 2013; 9(1):e1002866. https://doi.org/10.1371/journal.pcbi.1002866 PMID: 23382661

10. Szavits-Nossan J, Ciandrini L. Inferring efficiency of translation initiation and elongation from ribosome profiling. Nucleic acids research. 2020; 48(17):9478–9490. https://doi.org/10.1093/nar/gkaa678 PMID: 32821926

11. Tarantola A. Inverse Problem Theory and Methods for Model Parameter Estimation. USA: Society for Industrial and Applied Mathematics; 2004.

12. MacDonald CT, Gibbs JH, Pipkin AC. Kinetics of biopolymerization on nucleic acid templates. Biopolymers: Original Research on Biomolecules. 1968; 6(1):1–25. https://doi.org/10.1002/bip.1968.360060102 PMID: 5641411

13. Ibrahim Y, Dorignac J, Geniet F, Chevalier C, Walter JC, Walliser NO, et al. Chain decay and rates disorder in the totally asymmetric simple exclusion process. arXiv:2310.00329 [Preprint]. 2023 [posted 2023 Sep 30]. Available from: https://arxiv.org/abs/2310.00329
14. Heyer E, Moore M. Redefining the Translational Status of 80S Monosomes. Cell. 2016; 164(4):757-769. https://doi.org/10.1016/j.cell.2016.01.003 PMID: 26871635

15. Cox DR. Renewal theory. M. s. bartlett, f.r.s (methuen & co. ltd, science paperbacks) ed.; 1970.

16. D J Daley DVJ. An Introduction to the Theory of Point Processes. 2nd ed.; 2003. Available from: https://link.springer.com/book/10.1007/b97277.

17. Derrida B, Evans MR, Hakim V, Pasquier V. Exact solution of a 1D asymmetric exclusion model using a matrix formulation. Journal of Physics A: Mathematical and General. 1993; 26(7):1493-1517. https://doi.org/10.1088/0305-4470/26/7/011

18. Bernstein JA, Khodursky AB, Lin PH, Lin-Chao S, Cohen SN. Global analysis of mRNA decay and abundance in Escherichia coli at single-gene resolution using two-color fluorescent DNA microarrays. Proceedings of the National Academy of Sciences. 2002; 99(15):9697-9702. https://doi.org/10.1073/pnas.112318199 PMID: 12119387

19. Wang Y, Liu CL, Storey JD, Tibshirani RJ, Herschlag D, Brown PO. Precision and functional specificity in mRNA decay. Proceedings of the National Academy of Sciences. 2002; 99:5860-5865. https://doi.org/10.1073/pnas.092538799 PMID: 11972065

20. Sharova LV, Sharov AA, Nedorezov T, Piao Y, Shaik N, Ko MSH. Database for mRNA Half-Life of 19 977 Genes Obtained by DNA Microarray Analysis of Pluripotent and Differentiating Mouse Embryonic Stem Cells. DNA Research. 2009; 16(1):45-58. https://doi.org/10.1093/dnares/dsn030 PMID: 19001483

21. Schwanhäusser B, Busse D, Li N, Dittmar G, Schuchhardt J, Wolf J, et al. Global quantification of mammalian gene expression control. Nature. 2011; 473(7347):337-342. https://doi.org/10.1038/nature10098 PMID: 21593866

22. Young R, Bremer H. Polypeptide-chain-elongation rate in Escherichia coli B/r as a function of growth rate. The Biochemical Journal. 1976; 160(2):185-194. https://doi.org/10.1042/bj1600185 PMID: 795428

23. Karpins TV, Greenwood DJ, Sams CE, Ammons JT. RNA:protein ratio of the unicellular organism as a characteristic of phosphorous and nitrogen stoichiometry and of the cellular requirement of ribosomes for protein synthesis. BMC Biology. 2006; 4:30. https://doi.org/10.1186/1741-7007-4-30 PMID: 16953894

24. Ingolia N, Lareau L, Weissman J. Ribosome Profiling of Mouse Embryonic Stem Cells Reveals the Complexity and Dynamics of Mammalian Proteomes. Cell. 2011; 147(4):789-802. https://doi.org/10.1016/j.cell.2011.10.002 PMID: 22056041

25. Yan X, Hoek T, Vale R, Tanenbaum M. Dynamics of Translation of Single mRNA Molecules In Vivo. Cell. 2016; 165(4):976-989. https://doi.org/10.1016/j.cell.2016.04.034 PMID: 27153498

26. Trösemeier JH, Rudorf S, Loessner H, Hofner B, Reuter A, Schulenborg T, et al. Optimizing the dynamics of protein expression. Scientific Reports. 2019; 9(1). https://doi.org/10.1038/s41598-019-43857-5 PMID: 31101858

27. Relier S, Ripoll J, Guillorit H, Amalric A, Achour C, Boissière F, et al. FTO-mediated cytoplasmic m6A demethylation adjusts stem-like properties in colorectal cancer cell. Nature Communications. 2021; 12(1):1716. https://doi.org/10.1038/s41467-021-21758-4 PMID: 33741917

28. Gao X, Wan J, Qian SB. Genome-Wide Profiling of Alternative Translation Initiation Sites. Methods in Molecular Biology (Clifton, NJ). 2016; 1358:303-316. https://doi.org/10.1007/978-1-4939-3067-8_19 PMID: 26463392

29. Paulet D, David A, Rivals E. Ribo-seq enlightens codon usage bias. DNA Research. 2017; 24(3):303. https://doi.org/10.1093/dnares/dsw062 PMID: 28168289

30. Philippe N, Salson M, Commes T, Rivals E. CRAC: an integrated approach to the analysis of RNA-seq reads. Genome Biology. 2013; 14(3):R30. https://doi.org/10.1186/gb-2013-14-3-r30 PMID: 23537109

31. Stadler M, Fire A. Wobble base-pairing slows in vivo translation elongation in metazoans. RNA. 2011; 17(12):2063-2073. https://doi.org/10.1261/rna.02890211 PMID: 22045228

32. Heintz N, Sive HL, Roeder RG. Regulation of human histone gene expression: kinetics of accumulation and changes in the rate of synthesis and in the half-lives of individual histone mRNAs during the HeLa cell cycle. Molecular and Cellular Biology. 1983; 3(4):539-550. https://doi.org/10.1128/mcb.3.4.539 PMID: 6406835

33. scipy.optimize.minimize—SciPy v1.8.0 Manual; 2022. Available from: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html.

34. Ciandrini L, Neri I, Walter JC, Dauloudet O, Parmeggiani A. Motor protein traffic regulation by supply-demand balance of resources. Physical biology. 2014; 11(5):056006. https://doi.org/10.1088/1478-3975/11/5/056006 PMID: 25204752
35. Dauloudet O, Neri I, Walter JC, Dorignac J, Geniet F, Parmeggiani A. Modelling the effect of ribosome mobility on the rate of protein synthesis. The European Physical Journal E. 2021; 44(2):1-15. https://doi.org/10.1140/epje/s10189-021-00019-8

36. Dauloudet O, Neri I, Walter JC, Dorignac J, Geniet F, Parmeggiani A. Erratum to: Modelling the effect of ribosome mobility on the rate of protein synthesis. The European Physical Journal E. 2021; 44(4):1-1. https://doi.org/10.1140/epje/s10189-021-00061-6 PMID: 33900482

37. Hia F, Yang SF, Shichino Y, Yoshinaga M, Murakawa Y, Vandenbon A, et al. Codon bias confers stability to human mRNA s. EMBO reports. 2019; 20(11):e48220. https://doi.org/10.15252/embr.201948220 PMID: 31482640

38. Edri S, Tuller T. Quantifying the effect of ribosomal density on mRNA stability. PLoS One. 2014; 9(7): e102308. https://doi.org/10.1371/journal.pone.0102308 PMID: 25020060

39. Sin C, Chiarugi D, Valleriani A. Quantitative assessment of ribosome drop-off in E. coli. Nucleic acids research. 2016; 44(6):2528-2537. https://doi.org/10.1093/nar/gkw137 PMID: 26935582