Ribosome Flow Model with Extended Objects

TL;DR

The paper introduces the RFMEO, a deterministic model for ribosome flow along mRNA that incorporates realistic features and allows for rigorous analysis, improving understanding of translation dynamics.

Contribution

The RFMEO extends previous models by including ribosome size and non-homogeneous rates, and provides a mathematically tractable framework for analyzing translation.

Findings

RFMEO converges to a unique steady-state density.

Protein production rate reaches a steady-state.

Model entrains to periodic transition rates.

Abstract

We study a deterministic mechanistic model for the flow of ribosomes along the mRNA molecule, called the ribosome flow model with extended objects (RFMEO). This model encapsulates many realistic features of translation including non-homogeneous transition rates along the mRNA, the fact that every ribosome covers several codons, and the fact that ribosomes cannot overtake one another. The RFMEO is a mean-field approximation of an important model from statistical mechanics called the totally asymmetric simple exclusion process with extended objects (TASEPEO). We demonstrate that the RFMEO describes biophysical aspects of translation better than previous mean-field approximations, and that its predictions correlate well with those of TASEPEO. However, unlike TASEPEO, the RFMEO is amenable to rigorous analysis using tools from systems and control theory. We show that the ribosome density profile along the mRNA in the RFMEO converges to a unique steady-state density that depends on the length of the mRNA, the transition rates along it, and the number of codons covered by every ribosome, but not on the initial density of ribosomes along the mRNA. In particular, the protein production rate also converges to a unique steady-state. Furthermore, if the transition rates along the mRNA are periodic with a common period T then the ribosome density along the mRNA and the protein production rate converge to a unique periodic pattern with period T, that is, the model entrains to periodic excitations in the transition rates. We believe that the RFMEO could be useful for modeling, understanding, and re-engineering translation as well as other important biological processes.

Equations

$$\mu :={\left[\begin{array}{c}\alpha ,{\gamma }_1,\dots ,{\gamma }_{N-1},\beta \end{array}\right]}^{\prime }.$$

$${\dot{x}}_1$$

$${\dot{x}}_2$$

$${\dot{x}}_3$$

$$⋮$$

$${\dot{x}}_{n-1}$$

$${\dot{x}}_n$$

$${\dot{x}}_i=h_{i-1}(x)-h_i(x),i=1,\dots ,n,$$

$$C^n:=\{x\in {\mathbb{R}}^n:x_i\in [0,1],i=1,\dots ,n\}.$$

$$y_i(t):=\sum _{j=\max \{1,i-\ell +1\}}^i{x}_j(t),i=1,\dots ,n.$$

$$P=\left[\begin{array}{cccc}1& 0& 0& 0\\ 1& 1& 0& 0\\ 1& 1& 1& 0\\ 0& 1& 1& 1\end{array}\right].$$

$${\dot{x}}_i=q_{i-1}(x)-q_i(x),i=1,\dots ,n.$$

$$q_i(x):={\lambda }_i{x}_i(1-y_{i+\ell }),i=0,\dots ,n,$$

$${\dot{x}}_1$$

$$={\lambda }_0(1-y_{\ell })-{\lambda }_1{x}_1(1-y_{\ell +1}).$$

$$q_i(x)={\lambda }_i{x}_i(1-y_{i+3})={\lambda }_i{x}_i(1-x_{i+1}-x_{i+2}-x_{i+3});$$

$${\dot{x}}_j$$

$$={\lambda }_{j-1}{x}_{j-1}(1-y_{\ell +j-1})-{\lambda }_j{x}_j(1-y_{\ell +j})$$

$$={\lambda }_{j-1}{x}_{j-1}-{\lambda }_j{x}_j.$$

$${\dot{x}}_1$$

$$={\lambda }_0(1-x_1-x_2)-{\lambda }_1{x}_1(1-x_2-x_3),$$

$${\dot{x}}_2$$

$$={\lambda }_1{x}_1(1-x_2-x_3)-{\lambda }_2{x}_2(1-x_3-x_4),$$

$${\dot{x}}_3$$

$$={\lambda }_2{x}_2(1-x_3-x_4)-{\lambda }_3{x}_3,$$

$${\dot{x}}_4$$

$$x(0)={\left[\begin{array}{cccc}1& 0.1& 1& 1\end{array}\right]}^{\prime },$$

$${\dot{x}}_1(0)=0.1({\lambda }_1-{\lambda }_0).$$

$$J(x)=\left[\begin{array}{cccc}-{\lambda }_0-{\lambda }_1(1-x_2-x_3)& -{\lambda }_0+{\lambda }_1{x}_1& {\lambda }_1{x}_1& 0\\ {\lambda }_1(1-x_2-x_3)& -{\lambda }_1{x}_1-{\lambda }_2(1-x_3-x_4)& -{\lambda }_1{x}_1+{\lambda }_2{x}_2& {\lambda }_2{x}_2\\ 0& {\lambda }_2(1-x_3-x_4)& -{\lambda }_2{x}_2-{\lambda }_3& -{\lambda }_2{x}_2\\ 0& 0& {\lambda }_3& -{\lambda }_4\end{array}\right].$$

$${\dot{y}}_i=$$

$$-{\lambda }_i\left(\sum _{k=0}^{\lceil (i-\ell )/\ell \rceil }(y_{i-k\ell }-y_{i-k\ell -1})\right)(1-y_{i+\ell }),$$

$${\dot{y}}_1$$

$${\dot{y}}_2$$

$${\dot{y}}_3$$

$${\dot{y}}_4$$

$$H:=\{x\in {\mathbb{R}}^n:x\in C^n\text{ and }Px\in C^n\}.$$

$$x(t,a)\in H_{\tau },\text{ for all }t\ge \tau .$$

$$x_i(t,a),y_i(t,a)\in [\delta ,1-\delta ],\text{ for all }t\ge \tau ,$$

$$\dot{x}(t)=f(t,x(t)),$$

Full text

The RFMEO is a mean-field approximation of an important model from statistical mechanics called the totally asymmetric simple exclusion process with extended objects (TASEPEO). We demonstrate that the RFMEO describes biophysical aspects of translation better than previous mean-field approximations, and that its predictions correlate well with those of TASEPEO. However, unlike TASEPEO, the RFMEO is amenable to rigorous analysis using tools from systems and control theory. We show that the ribosome density profile along the mRNA in the RFMEO converges to a unique steady-state density that depends on the length of the mRNA, the transition rates along it, and the number of codons covered by every ribosome, but not on the initial density of ribosomes along the mRNA. In particular, the protein production rate also converges to a unique steady-state. Furthermore, if the transition rates along the mRNA are periodic with a common period $T$ then the ribosome density along the mRNA and the protein production rate converge to a unique periodic pattern with period $T$, that is, the model entrains to periodic excitations in the transition rates.

Analysis and simulations of the RFMEO demonstrate several counterintuitive results. For example, increasing the ribosome footprint may sometimes lead to an increase in the production rate. Also, for large values of the footprint the steady-state density along the mRNA may be quite complex (e.g. with quasi-periodic patterns) even for relatively simple (and non-periodic) transition rates along the mRNA. This implies that inferring the transition rates from the ribosome density may be non-trivial.

We believe that the RFMEO could be useful for modeling, understanding, and re-engineering translation as well as other important biological processes.

I Introduction

Gene expression is a multi-stage process for converting the information inscribed in the DNA to proteins. During the transcription stage, the information in the DNA of a specific gene is copied into messenger RNA (mRNA). In the translation stage, complex macro-molecules called ribosomes bind (at the initiation phase) to the mRNA and unidirectionally decode each codon (at the elongation phase) into the corresponding amino-acid that is delivered to the awaiting ribosome by transfer RNA (tRNA). Finally, at the termination phase, the ribosome detaches from the mRNA, the amino-acid sequence is released, folded, and becomes a functional protein (in some cases, post-translation modifications may occur) [1]. The output rate of ribosomes from the mRNA, which is also the rate in which proteins are generated, is called the protein translation rate, or production rate.

Translation occurs in all living organisms, and under almost all conditions. Thus, understanding the factors that affect translation has important implications to many scientific disciplines, including medicine, evolutionary biology, and synthetic biology. Deriving and analyzing mechanistic models of translation is important for developing a better understanding of this complex, dynamical, and tightly-regulated process. Such models can also aid in integrating and analyzing the rapidly increasing experimental findings related to translation (see, e.g., [10, 55, 54, 7, 47, 12, 40, 69]).

Mechanistic models of translation describe the dynamics of ribosome movement along the mRNA molecule, with parameters that encode the various translation factors affecting the codon decoding times along the mRNA molecule. Several such models have been suggested based on different paradigms ranging from Petri nets [5] to probabilistic Boolean networks [67]. For more details, see the survey papers [60, 69].

The totally asymmetric simple exclusion process (TASEP) [49, 68] is a fundamental model in non-equilibrium statistical mechanics that has been used to model numerous natural and artificial processes [46, 69], including ribosome flow during mRNA translation. In TASEP, particles stochastically hop between consecutive sites along an ordered lattice of $N$ sites. However, a particle cannot hop to an already occupied site. TASEP encapsulates both the unidirectional flow of ribosomes along the mRNA molecule, and the interaction between the particles, as a particle in site $i$ blocks the movement of a particle in site $i-1$. This hard exclusion principle models particles that have “volume” and thus cannot overtake one other. In the context of translation, the lattice represents the mRNA molecule, and the particles are the ribosomes. The rate of hoping from site $i$ to site $(i+1)$ is denoted by $\gamma_{i}$. A particle can hop to [from] the first [last] site of the lattice at a rate $\alpha$ [$\beta$]. The flow through the lattice converges to a steady-state value that depends on $N$ and the vector of parameters:

[TABLE]

The special case where all the internal hoping rates are assumed to be equal and normalized to one, i.e. $\gamma_{i}:=1$, $i=1,\dots,N-1$, is referred to as the homogeneous TASEP (HTASEP).

In TASEP a particle occupies a single site. However, in translation every ribosome occupies not only the codon it is translating, but also codons after and before it. More precisely, the ribosome footprint is about $10$ to $11$ codons, and its exit tunnel length is about $31$ codons [1, 23, 58, 66, 21]. In TASEP with extended objects (TASEPEO), a particle occupies multiple sites along the lattice [29, 30, 14, 51, 50, 27, 49]. For TASEPEO with open-boundary conditions (i.e. where the two sides of the lattice are connected to two particle reservoirs, as assumed here) few rigorous analytical results are known [14]. Mean-field approximations, domain-wall arguments, and extensive Monte Carlo simulations suggest that the homogeneous TASEPEO converges to a steady-state, and that the model has the same phase-diagram as HTASEP, i.e. it contains three phases: low-density, high-density, and maximal current. The phase boundaries depend on the extended object size [50]. TASEPEO with two types of object sizes was studied in [19]. It is important to mention that the extended objects concept is relevant for other intracellular processes e.g. transcription [16, 44, 26].

The ribosome flow model (RFM) [43] is a deterministic mathematical model for mRNA translation, obtained via a mean-field approximation of TASEP with open-boundary conditions. As such, it also inherits the property that the particle size is equal to the site size. When the RFM is used to model translation based on real biological data, this issue is handled by coarse-graining the mRNA molecule into sites composed of several consecutive codons. For every site the translation time of each codon in the site is used to determine the translation time of the site in the RFM (see e.g. [43]). It is not clear, however, how to systematically coarse-grain the mRNA in a way that yields the best fidelity between the model structure and parameters and the biological reality.

In this paper, we analyze for the first time a mean-field approximation of TASEPEO. This is a deterministic model that we refer to as the ribosome flow model with extended objects (RFMEO). Using the theory of contractive dynamical systems, we rigorously prove that the RFMEO always converges to a steady-state. In other words, the density profile of ribosomes along the mRNA molecule always converges to a unique steady-state, and thus so does the protein production rate. This shows that the RFMEO is robust in the sense that perturbations (e.g. due to stochastic noise in the biochemical reactions) in the ribosome density and production rate die out with time. This also means that we can reduce the problem of studying the density profile and protein production rate to studying the steady-state profile and production rate. We also prove that the RFMEO entrains (or frequency-locks) to periodic excitations. We show using simulations that the RFMEO, unlike the RFM, correlates well with TASEPEO.

The remainder of this paper is organized as follows. The next section briefly reviews the RFM. Section III describes the RFMEO. Section IV describes our main theoretical results on the properties of the RFMEO. Section V studies the correlation between RFMEO and TASEPEO. The final section summarizes and describes several directions for further research. To increase the readability of this paper, all the proofs are placed in Appendix A. Appendix B describes how the RFMEO can be derived by a mean-field approximation of TASEPEO.

II Ribosome Flow Model (RFM)

The RFM [43] is a deterministic model for mRNA translation that can be derived by a mean-field approximation of TASEP (see, e.g., [46, section 4.9.7], [4, p. R345] (see also Appendix B in the special case where the extended object size is equal to one site unit). In the RFM, mRNA molecules are coarse-grained into $n$ consecutive sites of codons. The state variable $x_{i}(t):\mathbb{R}_{+}\to[0,1]$, $i=1,\dots,n$, describes the normalized ribosomal occupancy level at site $i$ at time $t$, where $x_{i}(t)=1$ [$x_{i}(t)=0$] indicates that site $i$ is completely full [empty] at time $t$. The model includes $n+1$ positive parameters that describe the maximal possible transition rate between the sites: the initiation rate into the chain $\lambda_{0}$, the elongation (or transition) rate from site $i$ to site $(i+1)$ $\lambda_{i}$, $i=1,\dots,n-1$, and the exit rate $\lambda_{n}$.

The dynamics of the RFM with $n$ sites is given by $n$ nonlinear first-order ordinary differential equations:

[TABLE]

If we let $x_{0}(t):=1$ and $x_{n+1}(t):=0$, then (II) can be written more succinctly as

[TABLE]

where $h_{i}(x):=\lambda_{i}x_{i}(1-x_{i+1})$. This can be explained as follows. The flow of particles from site $i$ to site $(i+1)$ at time $t$ is $\lambda_{i}x_{i}(1-x_{i+1})$. This flow increases with the density at site $i$, and decreases as site $(i+1)$ becomes fuller. This corresponds to a “soft” version of a simple exclusion principle: since the particles have volume, the input rate to site $i$ decreases as the number of particles in that site increases. Note that the maximal possible flow from site $i$ to site $(i+1)$ is the transition rate $\lambda_{i}$. Thus Eq. (3) simply states that the change in the density at site $i$ at time $t$ is the input rate to site $i$ (from site $i-1$) minus the output rate (to site $i+1$) at time $t$.

The ribosome exit rate from site $n$ at time $t$ is equal to the protein production (or translation) rate at time $t$, and is denoted by $R(t):=\lambda_{n}x_{n}(t)$ (see Fig. 1). Note that $x_{i}$ is dimensionless, and that every rate $\lambda_{i}$ has units of 1/time.

Let $x(t,a)$ denote the solution of (II) at time $t\geq 0$ for the initial condition $x(0)=a$. Since the state-variables correspond to normalized occupancy levels, we always assume that $a$ belongs to the closed $n$-dimensional unit cube:

[TABLE]

Let $\operatorname{{\mathrm{i}nt}}(C^{n})$ denote the interior of $C^{n}$, and let $\partial C^{n}$ denote the boundary of $C^{n}$. It was shown in [34] that if $a\in C^{n}$ then $x(t,a)\in C^{n}$ for all $t\geq 0$, that is, $C^{n}$ is an invariant set of the dynamics. Ref. [34] also showed that the RFM is a tridiagonal cooperative dynamical system [52], and that this implies that (II) admits a unique steady-state point $e=e(\lambda_{0},\dots,\lambda_{n})\in\operatorname{{\mathrm{i}nt}}(C^{n})$, that is globally asymptotically stable, that is, $\lim_{t\to\infty}x(t,a)=e$, for all $a\in C^{n}$ (see also [31]). In particular, the production rate converges to the steady-state value $R:=\lambda_{n}e_{n}$.

An important advantage of the RFM (e.g. as compared to TASEP) is that it is amenable to mathematical analysis using various tools from systems and control theory. Furthermore, most of the analysis results hold for the general, nonhomogeneous case (i.e. when the transition rates all differ from one another). The RFM has been used to address many important biological problems including the sensitivity of the production rate to small changes in the transition rate, maximizing and minimizing the production rate in an optimal manner, analysis of the effect of competition for shared resources in translation, and more [33, 63, 34, 35, 31, 38, 39, 42, 62, 65, 41, 64].

To account for the fact that each ribosome covers several codons, we analyze here the RFMEO, which is a mean-field approximation of TASEPEO (see Appendix B for more details). An integer $\ell\geq 1$ describes the number of site units covered by each particle. The exclusion principle now implies that the rate of flow from site $i$ to site $(i+1)$ is $\lambda_{i}x_{i}(1-x_{i+1}-\cdots-x_{i+\ell})$. Indeed, since the particle covers the next $\ell$ sites, as the density in any of the $\ell$ consecutive sites increases the rate of movement slows down. Note that $\ell=1$ yields the RFM, so the RFM is a special case of the RFMEO.

Nevertheless, the RFMEO is a significant generalization of the RFM and its dynamics is quite different from that of the RFM. For example, the RFMEO, unlike the RFM, is not a cooperative system; it does not satisfy the particle-hole symmetry of the RFM (and of TASEP) [64, 4], and unlike the RFM, the RFMEO with $\ell>1$ is not a tridiagonal dynamical system.

III Ribosome Flow Model with Extended Objects (RFMEO)

Being a large complex of molecules, each ribosome typically covers between $10$ to $11$ codons and the geometry (e.g. length of the exit tunnel) can be longer than 30 codons [1]. A drawback of the RFM and other standard mean field models for translation is that, without additional processing such as coarse-graining, each ribosome (“particle”) is assumed to cover a single site.

The RFMEO allows modeling the flow of ribosomes where every ribosome covers $1\leq\ell\leq n$ site units. We assume, without loss of generality, that the ribosome is translating the left-most site it is covering, and refer to this part of the ribosome as the reader. A similar assumption is used in TASEPEO (see, for example, [14, 51, 50, 27, 49, 15]). Thus, the statement “the ribosome is at site $i$” means that: the reader is located at site $i$; the ribosome is translating site $i$; its corresponding transition rate is $\lambda_{i}$; and sites $i,\dots,i+\ell-1$, are covered by this ribosome. As we will show below, the dynamical equations describing the RFMEO (and thus all the theoretical results in this paper) are the same for any chosen reader location (e.g. choosing the reader at location $\ell/2$ results in exactly the same RFMEO equations).

Let $x_{i}(t)$ denote the (normalized) reader occupancy level at site $i$ at time $t$, and let $y_{i}(t)$ denote the (normalized) coverage occupancy level at site $i$ at time $t$, that is,

[TABLE]

Indeed, since every ribosome covers $\ell$ sites, any ribosome that is located up to $\ell$ sites left to site $i$ contributes to the total ribosome coverage at site $i$. The term “normalized” here means that each $x_{i}(t)$ and each $y_{i}(t)$ takes values in the interval $[0,1]$ for all $t\geq 0$. The value zero corresponds to completely empty, and one means completely full. We refer to $1-y_{i}(t)$ as the “space” or “vacancy” level at site $i$ at time $t$.

Note that (4) implies that $y(t)=Px(t)$, where $P$ is the lower triangular matrix with all entries zero, except for the entries on the main diagonal and $(\ell-1)$ diagonals below the main diagonal that are ones. For example, for $n=4$ and $\ell=3$:

[TABLE]

The dynamics of the RFMEO with $n$ sites is given by $n$ nonlinear first-order ordinary differential equations:

[TABLE]

Here $q_{i-1}$ is the flow into site $i$ and $q_{i}$ is the flow out of site $i$. The expression for this flow is given by

[TABLE]

with $x_{0}(t)\equiv 1$, and $y_{j}(t)\equiv 0$ for all $j>n$.

Eq. (7) implies that the reader flow from site $i$ to site $(i+1)$ is proportional to $\lambda_{i}$, to the occupancy levels of readers at site $i$, and to the “space” or “vacancy” level at site $i+\ell$ (see Fig. 2). In particular,

• •

As the number of readers at site $i$ increases, the flow from site $i$ increases. This follows the same reasoning as in the RFM.

• •

When a reader located at site $i$ moves to site $(i+1)$, the coverage occupancies at sites $i+1,i+2,\dots,i+\ell-1$ do not change. However, the reader’s tail end will now occupy a new site, which is site $(i+\ell)$.

• •

The “vacancy” level at site $(i+\ell)$ is $(1-y_{i+\ell})$, since $y_{i+\ell}$ denotes the total coverage at site $(i+\ell)$.

To explain (6), consider for example the equation for the change in the density at site $1$ given by

[TABLE]

The term $\lambda_{0}(1-y_{\ell})$ represents the entry rate into site $1$. Indeed, since the entering ribosome will cover sites $1,2,\dots,\ell$, this entry rate decreases with the coverage density $y_{\ell}=x_{1}+\dots+x_{\ell}$. (In the literature on TASEPEO this is referred to as the “complete-entry” flow [14]). The term $\lambda_{1}x_{1}(1-y_{\ell+1})$ is the flow from site $1$ to site $2$. This increases with the occupancy at site $1$ and, similarly, decreases with the coverage occupancy $y_{\ell+1}$.

Remark 1

As noted above, the assumption that the “reading head” is located at the left hand-side of the ribosome is arbitrary, but the RFMEO equations do not depend on this assumption. To demonstrate this, consider for example the case $\ell=3$ and the three possible locations for the reader: (1) Left-most site. In this case $y_{j}=x_{j-2}+x_{j-1}+x_{j}$, so

[TABLE]

(2) Middle site. In this case $y_{j}=x_{j-1}+x_{j}+x_{j+1}$, and $q_{i}(x)=\lambda_{i}x_{i}(1-y_{i+2})$, and this yields (8); (3) Right-most site. In this case $y_{j}=x_{j}+x_{j+1}+x_{j+2}$, and $q_{i}(x)=\lambda_{i}x_{i}(1-y_{i+1})$, again yielding (8). Thus (7) is invariant to the reader location.

Consider an index $j\geq n-\ell+2$. Then $\ell+j-1\geq n+1$, so

[TABLE]

Thus, the equation describing the flow in these last sites is a linear equation. The same phenomena takes place in TASEPEO, as a ribosome “reading” the last $\ell$ codons must be the last particle on the lattice, with no others to impede its progress. Therefore, it can move without hindrance toward the exit end. The exit rate in this context is referred to as the “incremental-exit” rate [14].

The output rate of ribosomes from the chain, which is the protein production (or translation) rate, is denoted by $R(t):=\lambda_{n}x_{n}$.

Note that in the special case $\ell=1$ we have $y_{i}=x_{i}$ for all $i=1,\dots,n$, and then (6) reduces to the RFM.

Example 1

Consider a RFMEO with dimension $n=4$ and particle size $\ell=2$. Then (6) yields

[TABLE]

In the RFM the set $C^{n}$ is an invariant set of the dynamics. This is no longer true for the RFMEO. For example for the initial condition

[TABLE]

Eq. (1) yields

[TABLE]

This implies that $x_{1}(0^{+})>1$ for $\lambda_{1}>\lambda_{0}$, so $x(0^{+})\not\in C^{n}$. $\square$

Remark 2

The Jacobian matrix of the dynamics (1) is

[TABLE]

Note that there are off-diagonal entries here whose sign may change with time (e.g. $-\lambda_{0}+\lambda_{1}x_{1}$). This implies that the RFMEO, unlike the RFM, is not a cooperative dynamical system.

It is useful to explicitly write the dynamics of the RFMEO in terms of the coverage state-variables (i.e. the $y$ state-vector). Recall that the proofs of all the results are placed in Appendix A.

Proposition 1

The coverage state-variables in the RFMEO satisfy:

[TABLE]

where $\lceil z\rceil$ denotes the smallest integer that is larger than or equal to $z$.

Example 2

Consider the RFMEO with $n=4$ sites and particle size $\ell=2$. Then (1) yields

[TABLE]

IV Theoretical Results

We begin by defining the relevant state-space for the RFMEO. If for some $i$ we have $y_{i+\ell}>1$ then $q_{i}(x):=\lambda_{i}x_{i}(1-y_{i+\ell})<0$. This represents a backward flow that according to current knowledge does not take place in ribosome movement. Thus, it is useful to define the state-space as the region where such a backward flow does not take place, i.e. both the $x_{i}$s and the $y_{i}$s are between zero and one. This leads to defining

[TABLE]

Note that $H$ is a compact and convex set.

Example 3

Consider the RFMEO with $n=3$ sites and particle size $\ell=2$. The sets $H$ and $C^{3}$ are depicted in Fig. 3. $\square$

Note also that for the RFMEO with $n=4$ and $\ell=2$ the initial condition $x(0)$ in (10) is not in $H$ as $y_{2}(0)=x_{1}(0)+x_{2}(0)>1$.

From here on we refer to any value $x\in H$ as a feasible value. This represents a state such that every reader density and every coverage density is between zero and one.

IV-A Invariance and persistence

The next result shows that the boundary of $H$, denoted $\partial H$, is “repelling” towards the interior of $H$. This means that if the dynamics is initiated with a feasible value that includes a reader/coverage density equal to the extremal value zero [one] then the dynamics will immediately change this to a value larger than zero [smaller than one].

Proposition 2

For any $a\in\partial H$ the solution of the RFMEO satisfies $x(t,a)\in\operatorname{{\mathrm{i}nt}}(H)$ for all $t>0$.

Note that this result implies in particular that $H$ is an invariant set of the dynamics. In other words, if all the reader and coverage densities are initiated with feasible values (i.e. values between zero and one) at time $t=0$ then they remain feasible for all time $t\geq 0$.

The next result shows that the solutions of the RFMEO are “persistent” in the sense that they enter and remain in a set that is uniformly separated from the boundary of $H$. Furthermore, this happens “immediately”.

Proposition 3

For any $\tau>0$ there exists a compact and convex set $H_{\tau}$ that is strictly contained in $H$ such that for any $a\in H$,

[TABLE]

Note that this implies that for any $\tau>0$ there exists $\delta=\delta(\tau)\in(0,1/2)$ such that

[TABLE]

for all $i$ and all $a\in H$. In other words, all the reader and coverage densities “immediately” become and remain uniformly separated from the extreme values zero and one. This is a technical property, but as we will see below it will be useful in the analysis of the asymptotic properties of the RFMEO.

IV-B Contraction

Contraction theory is a powerful tool for analyzing nonlinear dynamical systems. In a contractive system, trajectories that emanate from different initial conditions approach each other at an exponential rate, that is, the distance between any pair of trajectories, as a function of time, decreases at an exponential rate [28, 45, 2].

Consider the time-varying dynamical system

[TABLE]

whose trajectories evolve on a compact and convex set $\Omega\subset\mathbb{R}^{n}$.

For $t\geq t_{0}\geq 0$, and $a\in\Omega$, let $x(t,t_{0},a)$ denote the solution of (12) at time $t$ for the initial condition $x(t_{0})=a$. Recall that system (12) is said to be contartive on $\Omega$ w.r.t. a norm $|\cdot|:\mathbb{R}^{n}\to\mathbb{R}_{+}$ if there exists $\gamma>0$ such that

[TABLE]

for all $t_{2}\geq t_{1}\geq 0$ and all $a,b\in\Omega$. This means that any two trajectories approch each other at an exponential rate $\gamma$.

To apply contraction theory to the RFMEO, we require the following generalization of contraction with respect to a fixed norm that has been introduced in [32]. The time-varying system (12) is said to be contractive after a small overshoot and short transient (SOST) on $\Omega$ w.r.t. a norm $|\cdot|:\mathbb{R}^{n}\to\mathbb{R}_{+}$ if for each $\varepsilon>0$ and each $\tau>0$ there exists $\gamma=\gamma(\tau,\varepsilon)>0$ such that

[TABLE]

for all $t_{2}\geq t_{1}\geq 0$ and all $a,b\in\Omega$. Comparing this to (13), we see that here contraction “kicks in” after an arbitrarily small time transient $\tau$ and with an arbitrarily small overshoot $(1+\varepsilon)$.

The next result applies these ideas to the RFMEO. Let $|\cdot|_{1}:\mathbb{R}^{n}\to\mathbb{R}_{+}$ denote the $L_{1}$ norm, i.e. for $z\in\mathbb{R}^{n}$, $|z|_{1}=|z_{1}|+\dots+|z_{n}|$.

Proposition 4

The RFMEO is SOST on $H$ w.r.t. the $L_{1}$ norm, that is, for each $\varepsilon>0$ and each $\tau>0$ there exists $\gamma=\gamma(\tau,\varepsilon)>0$ such that

[TABLE]

for all $t\geq 0$ and all $a,b\in H$.

Roughly speaking, this means the following. Fix two initial feasible densities in the RFMEO and consider how the two corresponding densities along the mRNA evolve in time. Then these densities become “more similar” to each other at an exponential rate. In particular, the initial density is “quickly forgotten”.

Subsections IV-C and IV-E below describe important asymptotic properties of the RFMEO that follow from Prop. 4.

IV-C Global asymptotic stability

Write the RFMEO (6) as $\dot{x}=g(x)$, with $g:H\to\mathbb{R}^{n}$. Since the compact and convex set $H$ is an invariant set of this dynamical system, it contains a steady-state point $e=e(\lambda_{0},\dots,\lambda_{n})$. In other words, $g(e)=0_{n}$, where $0_{n}$ denotes a column vector of $n$ zeros, and $x(t,e)\equiv e$ for all $t\geq 0$. Prop. 2 implies that $e\in\operatorname{{\mathrm{i}nt}}(H)$. Using (14) with $b:=e$ yields the following result.

Corollary 1

The RFMEO admits a globally asymptotically stable steady-state point $e\in\operatorname{{\mathrm{i}nt}}(H)$, i.e.

[TABLE]

This means that trajectories corresponding to different initial conditions all converge to the unique steady-state point, that depends on the transition rates $\lambda_{i}$s, particle size $\ell$, and the length of the chain $n$, but not on the initial condition.

Example 4

Fig. 4 depicts the trajectories of an RFMEO with dimension $n=3$, particle size $\ell=2$, and rates $\lambda_{0}=1.0$, $\lambda_{1}=1.2$, $\lambda_{2}=0.8$, and $\lambda_{3}=0.4$, for four different initial conditions on the boundary of $H$. It may be seen that each trajectory immediately enters and remains in the interior of $H$, and converges to a unique steady-state point $e\in\operatorname{{\mathrm{i}nt}}(H)$. $\square$

The next example demonstrates the contraction property. Let $1_{n}$ denote the column vector of $n$ ones.

Example 5

Consider the RFMEO with dimension $n=7$, particle size $\ell=3$, and rates $\lambda_{i}=1-\frac{i}{50}$, $i=0,\dots,7$. In this case the unique steady-state point is (all numbers are to four digit accuracy):

[TABLE]

Fig. 5 depicts $r(t):=|x(t,a)-e|_{1}$ as a function of time for $t\in[0,70]$ for the initial condition $a=(3/20)1_{7}$. It may be seen that the $L_{1}$ distance between the trajectory and $e$ monotonically decreases to zero. It may also be seen that the rate of convergence varies with time. This makes sense because we can interpret the RFMEO as an RFM with time-varying transition rates (see the proof of Prop. 4 in Appendix A), and thus a time-varying contraction rate. $\square$

Corollary 1 implies that the coverage occupancy $y_{i}(t)$ at site $i$ converges to the unique steady-state value:

[TABLE]

Define the mean reader occupancy at time $t$ by $\rho(t):=\frac{1}{n}{\sum_{i=1}^{n}x_{i}(t)},$ and the mean coverage occupancy at time $t$ by $\rho^{c}(t):=\frac{1}{n}{\sum_{i=1}^{n}y_{i}(t)}.$ Note that this implies that $\lim_{n\to\infty}\rho^{c}(t)=\ell\rho(t)$. Then the mean reader occupancy converges to the unique steady-state value

[TABLE]

and the mean coverage occupancy converges to the unique steady-state value

[TABLE]

The next example demonstrates the contraction property using a S. Cerevisiae gene. Let $0_{n}$ denote the column vector of $n$ zeros.

Example 6

We consider the highly-expressed S. Cerevisiae gene YLR110C that encodes a cell wall mannoprotein, and contains $133$ codons (excluding the stop codon). We modeled it using a RFMEO with $n=133$ and $\ell=10$. The value $\lambda_{0}=1.33131$ [in units of mRNAs/sec] was estimated based on the ribosome density per mRNA levels, as this value is expected to be approximately proportional to the initiation rate when initiation is rate limiting [43, 33]. The elongation rates $\lambda_{1},\dots,\lambda_{n}$, were estimated using ribo-seq data for the codon decoding rates [11], normalized so that the median elongation rate of all S. cerevisiae mRNAs becomes $6.4$ codons per second [24]. The rates are depicted in the top panel of Fig. 6 as a function of $i$. To study the rate of contraction, we calculated $e$ in the RFMEO (shown in the middle panel of Fig. 6), and simulated the dynamical system to obtain $r(t)=|x(t,a)-e|_{1}$ with $a=0_{133}$, as a function of $t$ [in seconds]. Note that the initial condition $a=0_{133}$ represents an mRNA with no ribosomes. The bottom panel of Fig. 6 depicts the relative $L_{1}$ distance in percentage, that is,

[TABLE]

as a function of $t$ [in seconds]. In this case, $\rho=0.0534$. It may be seen that the relative distance is less than $20\%$ already after about $30$ seconds. We note that typical S. Cerevisiae mRNAs half-lives is in the order of tens of minutes (see, for example, [48, 61, 17]). This suggests that typically the ribosome density on S. Cerevisiae mRNAs is “very close” to its steady-state value.

IV-D Analysis of the steady-state

It is important to understand how the steady-state density $e$ depends on the parameters of the RFMEO. To study this we begin by deriving some equations for $e$. At steady-state (i.e. for $x=e$), the left-hand side of all the equations in (6) is zero (i.e. $\dot{x}_{i}=0$, $i=1,\dots,n$), so

[TABLE]

where $z_{j}:=0$, for all $j>n$. This yields (see (15))

[TABLE]

We can express these in terms of the (generally unknown value) $R$ as:

[TABLE]

and this yields

[TABLE]

and

[TABLE]

Example 7

Consider the RFMEO with dimension $n=6$ and particle size $\ell=3$. Then, the steady-state production rate $R$ satisfies

[TABLE]

where

[TABLE]

It is clear that solving (IV-D) is in general non-trivial. Nevertheless, it can be solved in closed-form in some very special cases. The next example demonstrates this.

Example 8

Consider a RFMEO with $n$ sites and with ribosome size $\ell=n$. Then (IV-D) becomes

[TABLE]

and this yields

[TABLE]

and

[TABLE]

where $z:=\sum_{0\leq i_{1}<i_{2}<\dots<i_{n}\leq n}\lambda_{i_{1}}\lambda_{i_{2}}\dots\lambda_{i_{n}}$. To understand this, assume in addition that $\lambda_{0}=\dots=\lambda_{n}=\lambda_{c}$, i.e. all the rates are equal with $\lambda_{c}$ denoting the common value. Then (22) and (23) yield $e_{i}=1/(n+1)$ for all $i$, and

[TABLE]

This means that when the ribosome size is equal to the chain size and all the rates are equal then the steady-state density at each site is identical. This makes sense, as every ribosome covers all the sites in the chain.

Another tractable case is when $\ell=n-1$ and $\lambda_{0}=\dots=\lambda_{n}=\lambda_{c}$. Then (IV-D) yields

[TABLE]

and this admits the solution

[TABLE]

and

[TABLE]

Note that here $e_{1}>e_{2}=e_{3}=\dots=e_{n}$. This makes sense, because if there is a ribosome with reader at site $\geq 2$ then the tail of this $(n-1)$-sites long ribosome is either at site $n$ or already out of the chain, and so there is no hindrance for its movement. $\square$

Eq. (IV-D) can also be used to prove theoretical results. The next result shows that increasing any of the $\lambda_{i}$s increases $R$. In other words, increasing any of the transition rates along the mRNA molecule increases the steady-state protein production rate.

Proposition 5

Consider the RFMEO with dimension $n$ and particle size $\ell$. Then $\frac{\partial}{\partial\lambda_{i}}R>0$, for $i=0,\dots,n$.

In the special case where all the rates are equal, i.e.

[TABLE]

where $\lambda_{q}$ denotes the common value, we refer to the RFMEO as the totally homogeneous RFMEO (THRFMEO). In this case, it is possible to say more about the steady-state occupancies.

Proposition 6

Consider the THRFMEO with dimension $n$ and particle size $\ell$. Then

[TABLE]

This means that the steady-state reader occupancies monotonically decrease between sites $1$ and $(n-\ell+1)$ and are equal at the last $\ell$ sites. This may partially explain the decrease in ribosome density observed along the coding sequences from the 5’ end to the 3’ end (see, for example, [9, 21]).

Example 9

The steady-state reader occupancy levels of the RFMEO with dimension $n=40$ are depicted in Fig. 7 for three particle sizes: $\ell=1$ (corresponding to the RFM), $\ell=2$, and $\ell=3$. It may be observed that the steady-state reader occupancies monotonically decrease along the chain until the last $\ell$ densities that are equal. The corresponding steady-state production rates are $R=0.2513$ for $\ell=1$; $R=0.1265$ for $\ell=2$; and $R=0.0851$ for $\ell=3$. $\square$

IV-D1 Effect of particle size

It is interesting to analyze how the steady-state occupancies and production rate depend on the particle size $\ell$. One might naturally expect the steady-state production rate in the RFMEO to decrease as the particle size $\ell$ increases. Indeed, roughly speaking one may think of increasing the particle size as replacing small cars traveling along a unidirectional traffic lane with large trucks thus leading to more congestion. This is demonstrated by the next example.

Example 10

The steady-state reader occupancy levels in the RFMEO with dimension $n=60$, and rates $\lambda_{0}=\cdots=\lambda_{40}=1$, and $\lambda_{41}=\cdots=\lambda_{60}=1/5$, for four different particle sizes: $\ell=1$ (i.e. the RFM), $\ell=2$, $\ell=4$, and $\ell=8$ are depicted in Fig. 8. Note that the steady-state occupancy levels decrease with $\ell$. The transition rates here decrease from the value $1$ to $1/5$ at site $40$. Thus, we except to see a “traffic jam” of ribosomes before this site. For $\ell=1$ this is indeed what happens. However, for $\ell>1$ much more complicated patterns evolve. The steady-state densities follow a complicated quasi-periodic behavior, with period $\ell$, even though there is no such periodicity in the rates. $\square$

Example 11

Fig. 9 depicts the steady-state production rate $R$, the steady-state mean reader occupancy $\rho$, and the steady-state mean coverage occupancy $\rho^{c}$ as a function of the particle size $\ell$, for a THRFMEO with dimension $n=100$ and $\lambda_{q}=1$. It can be observed that the steady-state production rate and the mean reader occupancy decrease with $\ell$, whereas the steady-state mean coverage occupancy increases with $\ell$.

It is interesting to compare these results to the homogeneous TASEPEO. In the thermodynamical limit (i.e. as $N\to\infty$), the homogeneous TASEPEO with particle size $\ell$, and with $\alpha=\beta=1$ is in the maximal current phase, where the steady-state output rate is $J=1/(1+\sqrt{\ell})^{2}$, the mean reader density is $1/(\sqrt{\ell}(1+\sqrt{\ell}))$, and the mean coverage density is $\sqrt{\ell}/(1+\sqrt{\ell)}$ [50, 14]. Note that this implies that as $\ell$ goes to infinity the current and the mean reader density go to zero, whereas the mean coverage density goes to one. This is consistent with the results for the THRFMEO depicted in Fig. 9.

Fig. 10 depicts the steady-state production rate $R$ as a function of the particle size $\ell$ for a RFMEO with $n=100$, $\lambda_{0}=0.1$, and $\lambda_{i}=1$, $i=1,\dots,100$. In this case $\lambda_{0}$ is the rate limiting factor, and thus less “traffic jams” occur relative to the case $\lambda_{0}=1$. It may be seen that also in this case $R$ monotonically decreases with $\ell$. $\square$

The next result shows that for fixed rates the steady-state production rate in the RFMEO with $\ell>1$ is always smaller than the steady-state production rate in the RFMEO with $\ell=1$ (i.e. the RFM).

Proposition 7

Consider an RFMEO with dimension $n$, particle size $\ell>1$, and rates $\lambda_{i}$, $i=0,\dots,n$, admitting a steady-state production rate $R$. Consider also an RFM with the same dimension $n$, and the same rates $\lambda_{i}$, $i=0,\dots,n$, admitting a steady-state production rate $\bar{R}$. Then $R<\bar{R}$.

In many organisms longer genes have lower protein levels [10, 18]. There are many explanations and variables that may contribute to this correlation. However, is it possible that the relations between particle size and translation rate may have a (small) contribution to this correlation? It is possible that longer coding regions are related to longer proteins emerging from the ribosome during translation thus practically increasing the effective ribosome size. This hypothesis may be studied in synthetic system in the future.

Surprisingly, however, increasing $\ell$ does not always lead to a reduction in the production rate.

Example 12

Consider an RFMEO with dimension $n=3$, and rates

[TABLE]

We consider two cases $\ell=2$ and $\ell=3$, and for the sake of clarity we denote the steady-state values in the latter case by overbars. For $\ell=2$,

[TABLE]

For $\ell=3$,

[TABLE]

Thus, here increasing $\ell$ from $2$ to $3$ increases the production rate. To explain this, recall that in general increasing $\ell$ decreases the steady-state reader densities (see Figs. 7 and 8). This is also what happens here. Indeed,

[TABLE]

At steady-state, the entry rate into the chain is equal to the production rate, so $R=\lambda_{0}(1-e_{1}-e_{2})$ and $\bar{R}=\lambda_{0}(1-\bar{e}_{1}-\bar{e}_{2}-\bar{e}_{3})$. Since this is proportional to one minus the sum of densities, $\bar{R}>R$. Thus, in this particular case the increase in $\ell$ yields an increase in the production rate.

Similarly, it follows from (24) and (26) that for any $n>3$ increasing $\ell$ from $n-1$ to $n$ in the THRFMEO leads to an increase in the steady-state production rate. $\square$

IV-E Entrainment

Assume now that some or all of the transition rates $\lambda_{i}$ are not constants, but time-varying periodic functions with a common period $T$. In the context of translation, this corresponds for example to the case where the tRNA abundances vary in a periodic manner, with a common period $T$. More precisely, we say that a function $f$ is $T$-periodic if $f(t+T)=f(t)$ for all $t$. Assume that the transition rates are time-varying functions satisfying:

There exist $\delta_{1},\delta_{2}\in\mathbb{R}_{++}$ such that $0<\delta_{1}\leq\lambda_{i}(t)\leq\delta_{2}$, for all $i=0,\dots,n$, and all $t\geq 0$. 2. 2.

There exists a minimal $T>0$ such that every $\lambda_{i}(t)$ is a $T$-periodic function.

We refer to this case as the periodic RFMEO (PRFMEO). Note that the PRFMEO includes in particular the case where some of the rates are constant, as a constant function is $T$-periodic for every $T$. However, condition $2)$ above implies that the case where all the rates are constant is ruled out, as then the minimal $T$ is zero. Indeed, this case is just the RFMEO analyzed above.

The next result follows from combining the fact that the RFMEO is SOST on $H$ with known results on the entrainment of contractive systems to a periodic excitation (see e.g. [45]).

Theorem 1

The PRFMEO admits a unique function $\phi(\cdot):\mathbb{R}_{+}\to\operatorname{{\mathrm{i}nt}}(H)$, that is $T$-periodic, and for any $a\in H$ the trajectory $x(t,a)$ converges to $\phi$ as $t\to\infty$.

In other words, the PRFMEO admits a unique periodic solution, with period $T$, and every trajectory of the PRFMEO converges to this periodic solution. This means that the densities along the mRNA, and thus also the production rate entrain to the periodic excitation induced by the transition rates.

As a side note, we point that the RFMEO can also be used to model vehicular traffic. If traffic lights that change periodically produce the transition rates then Thm. 1 implies that the traffic density converges to a periodic pattern with the same period, i.e. the “green wave” concept (see, e.g., [25]).

The next example illustrates the dynamical behavior of the PRFMEO.

Example 13

Consider an PRFMEO with dimension $n=4$, particle size $\ell=2$, and transition rates

[TABLE]

Note that all the rates here are periodic, with a minimal common period $T=8$. Figure 11 depicts $x_{i}(t)$, $i=1,\dots,4$, as a function of $t$ for the initial condition $x(0)=\begin{bmatrix}0.2&0.2&0.2&0.2\end{bmatrix}^{\prime}$. It may be seen that each state-variable converges to a periodic function with period $T=8$. $\square$

IV-F Rate limiting steps in the RFM and the RFMEO

It has been shown that depending on the biological conditions and the specific organism both initiation and elongation may be rate limiting [56, 8, 69, 40, 53, 60, 22]. Since the RFMEO is a better model for biological translation than the RFM, it is interesting to study the rate limiting step in these two models. We now show that the transition from the low density phase, when initiation is rate limiting, to the high density phase, when elongation is rate limiting is different in the two models: in the RFMEO this transition will take place for a lower initiation rate.

We modeled four S. cerevisiae genes: YMR123W, YNL303W, YJR094W-A, and YBL094C using both an RFMEO with $\ell=10$ and an RFM, and considered the steady-state production rate and the steady-state mean density as a function of the initiation rate $\lambda_{0}$.

As was done in Example 6, the elongation rate $\lambda_{i}$ at each site, for both the RFMEO and the RFM, was estimated using ribo-seq data for the codon decoding rates normalized so that the median elongation rate of all S. cerevisiae mRNAs becomes $6.4$ codons per second. The site rate is simply the corresponding codon rate. These rates thus depend on various factors including availability of tRNA molecules, amino acids, Aminoacyl tRNA synthetase activity and concentration, and local mRNA folding [11, 1, 56].

Fig. 12 depicts the steady-state production rate as a function of $\lambda_{0}$ for the four S. cerevisiae genes for both the RFMEO (upper figure) and the RFM (lower figure). It may be seen that the transition from an initiation rate limiting stage to an elongation rate limiting stage occurs for a lower initiation value in the RFMEO as compared to the RFM.

Fig. 13 depicts the steady-state mean density as a function of $\lambda_{0}$ for the four genes and two models. Again, it can be seen that the transition from an initiation rate limiting stage to the elongation rate limiting stage occurs at lower initiations value in the RFMEO as compared to the RFM. This holds for all four genes.

V High correlation between RFMEO and TASEPEO

In this section, we show that the RFMEO correlates better with TASEPEO than the RFM, supporting the modeling of intracellular process with multi-site biological machines such as translation and transcription using the RFMEO.

The simulations of TASEPEO with dimension $N$, rates $\mu$ (see (1)), and particle size $\ell$ use a parallel update mode. At each time tick $t_{k}$, the sites along the lattice are scanned from site $N$ backwards to site $1$. If it is time to hop, and the site that is $\ell$ sites in front is empty then the reader advances to the consecutive site. If the site that is $\ell$ sites in front is occupied, the next hopping time, $t_{k}+\varepsilon_{k}$, is generated randomly. For site $i$, $\varepsilon_{k}$ is exponentially distributed with parameter $(1/\mu_{i+1})$ (see (1)). The occupancy at each site is averaged throughout the simulation, with the first $700,000$ cycles discarded in order to obtain the steady-state value. We use $\varrho\in\mathbb{R}^{N}_{+}$ to denote the steady-state reader density, $J:=\beta\varrho_{N}$ to denote the steady-state current (or output rate), and $\sigma:=(1/N)\sum_{i=1}^{N}\varrho_{i}$ for the steady-state mean reader density.

In the examples below, we numerically calculated the Pearson correlation coefficients between the steady-states of the RFMEO, TASEPEO, and RFM.

Example 14

Consider the RFMEO with dimension $n=75$, and transition rates $\lambda_{0}=\cdots=\lambda_{75}=1$. Let $\tilde{e}$ denote the steady-state density of an RFM with the same dimension and rates. We also simulated TASEPEO with dimension $N=75$ and rates $\mu=\lambda$. Fig. 14 depicts the Pearson correlation coefficient $r(e,\varrho)$ between the steady-state reader densities of the RFMEO and the TASEPEO, and the Pearson correlation coefficient $r(\tilde{e},\varrho)$ between the steady-state reader densities of the RFM and TASEPEO, as a function of $\ell\in\{1,\dots,30\}$. The corresponding p-values were all less than $10^{-50}$. It may be seen that $r(e,\varrho)$ and $r(\tilde{e},\varrho)$ are somewhat similar for $\ell\in\{1,\dots,5\}$, however for all $\ell>5$, $r(e,\varrho)>0.94$ whereas $r(\tilde{e},\varrho)$ decreases with $\ell$, and is equal to about $0.825$ for $\ell=30$. Of course, this makes sense as the RFMEO is a mean field approximation of TASEPO. $\square$

The following examples consider the case of non-homogeneous transition rates.

Example 15

Consider the RFMEO with dimension $n=40$, particle size $\ell=15$, and rates $\lambda_{i}=1+0.3\sin(2\pi i/41)$, $i=0,\dots,40$. Fig. 15 depicts the RFMEO steady-state reader density $e$, the TASEPEO steady-state reader density $\varrho$ for $\mu=\lambda$, and particle size $15$, and the steady-state density $\tilde{e}$ in the RFM with the same dimension and rates. It can be seen that $e$ provides a far better estimate of $\varrho$ than $\tilde{e}$. $\square$

In order to verify that the high correlation between RFMEO and TASEPO holds for a large set of parameters, we also simulated the case where the rates are drawn randomly.

Example 16

Consider the RFMEO with dimension $n=100$, particle size $\ell=10$, and rates

[TABLE]

where $\theta_{i}\sim U[-1/2,1/2]$ is a random variable uniformly distributed in the interval $[-1/2,1/2]$. We compared the steady-state production rates of this RFMEO with those of the corresponding TASEPO, and with two RFMs. One RFM with the same dimension and rates. Another RFM, that we refer to as RFM10, is an approximation of the chain with $10$ “codons”/site. Thus, it has dimension $(100/10)-1=9$, where each site contains $10$ consecutive sites of the RFMEO (other than the last site which contains the last $11$ consecutive sites of the RFMEO). The rates of RFM10 are $\gamma_{i}=(\sum_{k=10i}^{T_{i}}\lambda_{k}^{-1})^{-1}$, where $T_{i}=(10(i+1)-1)$ if $i<9$, and otherwise $T_{i}=100$. Note that since the dimension of this RFM$10$ is nine, it cannot be used to estimate the entire density profile of the TASEPEO with dimension $100$.

We ran $300$ tests, where in each test a new set of rates were drawn according to (29). Fig. 16 depicts the correlation between the steady-state production rates of (1) RFMEO and TASEPEO; (2) RFM (i.e. RFM with dimension $100$ and rates $\lambda_{i}$) and TASEPEO, and (3) RFM$10$ and TASEPEO, over the $300$ tests. It may be seen that the RFMEO provides the best correlation with TASEPEO.

Fig. 17 depicts the correlations between the steady-state mean densities for the same three cases. It may be seen that again the correlation between the RFMEO and TASEPEO is high ($r\simeq 0.927$). The correlation between the RFM$10$ and TASEPEO is slightly better ($r\simeq 0.944$), however, as stated above, RFM10 cannot be used to provide an estimate to the actual (per codon) density profile. $\square$

VI Discussion

We studied a deterministic mechanistic model for mRNA translation, the RFMEO, that encapsulates many realistic features of this biological process including the fact that every ribosome covers several codons and that ribosomes cannot overtake one another.

The RFMEO is a mean-field approximation of TASEPEO (see Appendix B) and, as demonstrated above, its simulation results often correlate well with those of TASEPEO. However, unlike TASEPEO, the RFMEO is amenable to rigorous analysis using tools from systems and control theory.

We proved that the RFMEO converges to a unique state-state density and steady-state production rate for any set of feasible transition rates. We follow the terminology used in physics, where an equilibrium point [steady-state] is characterized by a zero [constant but nonzero] total flow of energy [6]. The convergence to this unique steady-state takes place at an exponential rate. In this respect, the RFMEO is robust to the initial conditions.

One may naturally ask whether biological systems are at steady-state (that maybe more general than the steady-state here, e.g. a periodic trajectory). Models with a steady-state (or several of steady-states) have been found to be useful in numerous studies in systems biology (see, e.g. [20] and the references therein). In practice the state of the art routine biological experiments and their interpretation assume steady state as they are performed in a very specific experimental environment which is kept constant during the entire experiment (see, for example, [37, 59, 3]).

In particular, the steady-state in the RFM has been used to accurately predict several features of gene expression (see, e.g., [43, 69, 17]). Here, we used the RFMEO to model a highly-expressed S. cerevisiae gene. The rates were estimated based on biological data. In the resulting RFMEO the convergence to a state close to the steady-state takes approximately 30 seconds, whereas the mRNA half-life is of the order of tens of minutes. This suggests that at least in this case the steady-state assumption is justified.

An important question is how does the steady-state depend on the RFMEO parameters. We proved that increasing any of the RFMEO rates can only increase the steady-state production rate, and that in the totally homogeneous case (i.e. when all the rates are equal) the reader ribosomal density monotonically decreases along the mRNA. In addition, we proved that if one or more of the RFMEO rates are time-varying periodic functions, with a common period $T$, then the densities along the mRNA, and thus also the production rate converge to a periodic solution with period $T$.

The results reported here can shed light on various biophysical aspects of translation, and can be further studied experimentally. For example, our analysis suggests that higher decoding rates at the last $\ell$ codons of the coding region can be expected (since in this region no downstream ribosome can block the ribosome movement). This can be validated experimentally for example based on approaches that track the movement of ribosomes at high resolution [57].

In addition, analysis and simulations of the RFMEO demonstrate several surprising and counterintuitive results. For example, increasing the particle size $\ell$ (i.e. the ribosome footprint) may some times lead to an increase in the production rate. Also, for large $\ell$ the steady-state density along the mRNA may be quite complex (e.g. with quasi-periodic patterns) even for relatively simple (and non-periodic) transition rates. It will be interesting to see if similar patterns are observed experimentally by possibly engineering the codon elongation rates of heterologous or endogenous genes and monitoring translation [57, 21].

We believe that the RFMEO could be useful for modeling, understanding, and re-engineering translation. Specifically, the advantages of the model mentioned above should make it a better candidate than other alternative models for solving some of the open questions in the field [69].

An important topic for future research is using the RFMEO to model ribosome flow based on biological data. This is a challenging task, as many aspects of translation are still not clear. For example, translation initiation is affected by complex phenomena such as the number of free ribosomes, mRNA folding near the 5’end of the mRNA, UTR length and other features, the nucleotide composition surrounding the start codon, and more. In addition, current techniques for measuring ribosome densities provide partial, noisy, and biased data (see, for example, [13]). Thus, using them to estimate the parameters in a computational model like the RFMEO is a non trivial challenge.

Another research topic is using the RFMEO (and networks of RFMEOs) to study various phenomena such as competition for resources in mRNA translation [41, 69], transcription [16], and evolution of transcripts [69].

Acknowledgments

We are grateful to the anonymous referees for their comments that greatly helped in improving this paper.

Appendix A: Proofs

Proof of Prop. 1. Combining (4) and (6) yields

[TABLE]

By the definition of $y_{i}$, $x_{i}=y_{i}-y_{i-1}+x_{i-\ell},$ and iterating this yields

[TABLE]

Substituting this in (Appendix A: Proofs) yields (1). ∎

Proof of Prop. 2. Consider the RFMEO with $x(0)\in\partial H$. Then $y(0)=Px(0)$, and there exists an index $i$ such that either $x_{i}(0)\in\{0,1\}$ or $y_{i}(0)\in\{0,1\}$ and all the other entries of $x(0)$ and $y(0)$ are between zero and one. The proof is based on computing the derivatives of the state-variables at time zero, and showing that state-variables that are zero [one] become strictly larger than zero [strictly smaller than one] at time $0^{+}$. We assume throughout that $\ell\geq 2$, as otherwise the RFMEO reduces to the RFM and then the proof follows from the results in [31]. We consider several cases.

Case 1. Suppose that $y_{\ell}(0)=0$. This implies in particular that $x_{\ell}(0)=0$. By (Appendix A: Proofs),

[TABLE]

Thus, $y_{\ell}(0^{+})>0$. Note that this calculation also implies that for any $\tau>0$ there exists $\varepsilon_{\ell}=\varepsilon_{\ell}(\tau)>0$ such that $y_{\ell}(t,a)\geq\varepsilon_{\ell}$ for all $t\geq\tau$ and all $a\in H$.

Case 2. Suppose that $y_{\ell+1}(0)=0$. This implies in particular that $x_{\ell+1}(0)=0$, so $y_{\ell}(0)=y_{\ell}(0)-y_{\ell+1}(0)=x_{1}(0)-x_{\ell+1}(0)=x_{1}(0)$. By (Appendix A: Proofs),

[TABLE]

Combining this with the result in Case 1 implies that for any $\tau>0$ there exists $\varepsilon_{\ell+1}=\varepsilon_{\ell+1}(\tau)>0$ such that $y_{\ell+1}(t,a)\geq\varepsilon_{\ell+1}$ for all $t\geq\tau$ and all $a\in H$.

Continuing in this fashion shows that for any $\tau>0$ there exists $\varepsilon=\varepsilon(\tau)>0$ such that $y_{i}(t,a)\geq\varepsilon$ for all $i\in\{\ell,\ell+1,\dots,n\}$, all $t\geq\tau$, and all $a\in H$.

Case 3. Suppose that $x_{j}(0)=0$ for some $j$. Then there exists a minimal index $i$ such that $x_{i}(0)=0$. If $i=n$ then (6) yields

[TABLE]

By the definition of $i$, $x_{n-1}(0)>0$ and thus $x_{n}(0^{+})>0$.

Now suppose that $i=n-1$. Then (6) yields

[TABLE]

By the definition of $i$, $x_{n-2}(0)>0$. If $\ell>2$ then $1-y_{n+\ell-2}(0)=1$, and thus $x_{n-1}(0^{+})>0$. If $\ell\leq 2$ then $1-y_{n+\ell-2}(0)=1-y_{n}(0)=1-x_{n}(0)-x_{n-1}(0)=1-x_{n}(0)$. Thus, if $x_{n}(0)<1$ then $x_{n-1}(0^{+})>0$. Consider the case $x_{n}(0)=1$. Then $y_{n}(0)=x_{n-1}(0)+x_{n}(0)=1$, so

[TABLE]

This means that $y_{n}(0^{+})<1$, so again we conclude that $x_{n-1}(0^{+})>0$.

Continuing in this fashion shows that if $x_{j}(0)=0$ for some $j$ then $x_{j}(0^{+})>0$. The analysis in all the other relevant cases is very similar, and thus omitted. ∎

Proof of Prop. 3. This follows from the fact that $H$ is compact, convex and with a repelling boundary; see [36, Thm. 2] (see also [32]). ∎

Proof of Prop. 4. Pick $\varepsilon,\tau>0$ and $a,b\in H$. By Prop. 3, there exists $\delta=\delta(\tau)\in(0,1/2)$ such that for all $i$ and all $t\geq\tau$,

[TABLE]

Write the $q_{j}$s in (6) as

[TABLE]

where $\eta_{j}(t):=\lambda_{j}\frac{1-y_{j+\ell}(t)}{1-x_{j+1}(t)}.$ Note that (32) implies that

[TABLE]

for all $j$ and all $t\geq\tau$. Using this notation, the RFMEO in (6) can be written as the time-varying system

[TABLE]

This means that for all $t\geq\tau$ the RFMEO can be interpreted as an RFM with time-varying transition rates $\eta_{j}(t)$ that, by (33), are uniformly bounded and uniformly separated from zero for all $t\geq\tau$. Now the results in [31] imply that there exists $\gamma:=\gamma(\varepsilon)$ such that after time $\tau$ the solutions are contractive with overshoot $(1+\varepsilon)$, and this completes the proof. ∎

Proof of Prop. 5. Consider two RFMEOs, both with the same dimension $n$ and particle size $\ell$. The first with rates $\lambda_{0},\dots,\lambda_{n}$, admits a steady-state density $e$, and a steady-state production rate $R$, and the second with rates $\tilde{\lambda}_{0},\dots,\tilde{\lambda}_{n}$, admits a steady-state density $\tilde{e}$ and a steady-state production rate $\tilde{R}$. Assume that there exists an index $j\in\{0,\dots,n\}$, such that $\tilde{\lambda}_{i}=\lambda_{i}$ for all $i\neq j$, and

[TABLE]

We need to show that $\tilde{R}>R$. Seeking a contradiction, assume that

[TABLE]

We start with the case $j=n$. Combining (35), (34) and (IV-D) implies that $\tilde{e}_{n}<e_{n}$, and $\tilde{e}_{n-k}\leq e_{n-k}$, $k=1,\dots,\ell-1$. This means that $\tilde{y}_{n}<y_{n}$, and combining this with (35) and (IV-D) implies that $\tilde{e}_{n-\ell}<e_{n-\ell}$, and so $\tilde{y}_{n-1}<y_{n-1}$. Continuing in this way yields $\tilde{e}_{j}<e_{j}$, $j=1,\dots,n-\ell$. In particular, $\tilde{e}_{1}+\dots+\tilde{e}_{\ell}<e_{1}+\dots+e_{\ell}$, and using (IV-D) results in $\tilde{R}>R$. This contradicts (35), and so we conclude that $\tilde{R}>R$ in the case where $\tilde{\lambda}_{n}>\lambda_{n}$.

Using the same approach for any $j\in\{0,\dots,n\}$, while combining the assumption in (35) with (34) and (IV-D), yields

[TABLE]

If $j>0$ then using $k=\ell$ in (Appendix A: Proofs) yields $\tilde{e}_{1}+\dots+\tilde{e}_{\ell}<e_{1}+\dots+e_{\ell}$, thus $\tilde{R}>R$, contradicting (35). If $j=0$ then using $k=\ell$ in (Appendix A: Proofs) yields $\tilde{e}_{1}+\dots+\tilde{e}_{\ell}\leq e_{1}+\dots+e_{\ell}$, but since $\tilde{\lambda}_{0}>\lambda_{0}$, this again yields $\tilde{R}>R$, contradicting (35). We conclude that $\tilde{R}>R$. ∎

Proof of Prop. 6. Consider (IV-D) with $\lambda_{0}=\cdots=\lambda_{n}$. Then

[TABLE]

Since $e_{n-\ell}(1-e_{n-\ell+1}-\cdots-e_{n})=e_{n-\ell}(1-z_{n})=e_{n-\ell+1}$, and $z_{n}\in(0,1)$, it follows that

[TABLE]

and combining this with (37) implies that

[TABLE]

Now, since $e_{n-\ell-1}(1-e_{n-\ell}-\cdots-e_{n-1})=e_{n-\ell-1}(1-z_{n-1})=e_{n-\ell}(1-z_{n})$, using (39) and the fact that $z_{n-1},z_{n}\in(0,1)$ imply that $e_{n-\ell-1}>e_{n-\ell}$ and thus $z_{n-2}>z_{n-1}$. Continuing in this way completes the proof. ∎

Proof of Prop. 7. Let $e$ $[\bar{e}]$ denote the steady-state reader density in the RFMEO [RFM]. We need to show that $\bar{R}>R$. Seeking a contradiction, assume that

[TABLE]

Combining this with (IV-D) for both the RFMEO with particle size $\ell$ and with particle size one (i.e. the RFM), it follows that $\lambda_{0}(1-\bar{e}_{1})\leq\lambda_{0}(1-e_{1}-\dots-e_{\ell})$, thus

[TABLE]

and since $e\in\operatorname{{\mathrm{i}nt}}(H)$ this yields

[TABLE]

Using (IV-D), (40), and (41), it follows that

[TABLE]

and since $e\in\operatorname{{\mathrm{i}nt}}(H)$ this yields

[TABLE]

Continuing in this way yields

[TABLE]

so in particular,

[TABLE]

On the other-hand using (40) and comparing the last $\ell$ equations in (IV-D) for both the RFMEO with particle size $\ell$ and with particle size one (i.e. the RFM), yields

[TABLE]

Now combining (Appendix A: Proofs) with (43) yields

[TABLE]

However, since $\bar{e}\in\operatorname{{\mathrm{i}nt}}(H)$, the term on the left-hand side here must be strictly positive. This contradiction completes the proof. ∎

Appendix B: RFMEO as a Mean-Field Approximation of TASEPEO

In this appendix, we show how the RFMEO can be derived from TASEPEO. We use a notation that is standard in the TASEPEO literature.

Consider TASEPEO with $N$ sites, rates $\mu$ defined in (1), extended object size $\ell$, and under the assumption that the reader is located at the left-most site of the object. Following MacDonald et. al. [29] (see also [30]) the current from site $i$ to site $i+1$ at time $t$ is given by (for simplicity we ignore boundary cases):111Note that in [29] the reader is defined to be in the right-most site of the object, and thus there the current is proportional to the probability that site $i$ has a reader and site $i+1$ is empty.

[TABLE]

where $\operatorname{Pr}(a)$ [$\operatorname{Pr}(a|b)$] denotes the probability of event $a$ [the conditional probability of event $a$ given event $b$] at time $t$. Since the conditional probability in (Appendix B: RFMEO as a Mean-Field Approximation of TASEPEO) is difficult to estimate, we apply what [14] calls a naive mean-field approximation, and replace (Appendix B: RFMEO as a Mean-Field Approximation of TASEPEO) by:

[TABLE]

We approximate the probabilities above by averaging the binary reader occupancies over an ensemble of TASEPEO systems, i.e. we replace $\operatorname{Pr}(\text{site }i\text{ has a reader})$ by $\rho^{r}_{i}(t):=\langle r_{i}(t)\rangle$, where $r_{i}(t)\in\{0,1\}$ is the reader occupancy at site $i$ at time $t$, and the operator $\langle\rangle$ denotes an average over the ensemble. This yields

[TABLE]

The change in the average reader occupancy at site $i$ at time $t$ is given by [29]:

[TABLE]

Introducing the notation $x_{i}(t):=\rho^{r}_{i}(t)$ and $\lambda_{i}:=\gamma_{i},$ we see that $J_{i\to i+1}(t)$ corresponds to $q_{i}(x)$ in (7), and (49) corresponds to (6) (see (4)). Thus, we obtained the RFMEO. In particular, the case $\ell=1$ in (48) corresponds to the dynamical equations of the RFM (see (3)).

At steady-state, we expect every $\rho^{r}_{i}(t)$ in TASEPEO to converge to, say, $\rho^{r}_{i}$, and then the currents between any two consecutive sites are all equal (but we are not aware of any rigorous proof of convergence in TASEPEO). The derivation above (including the boundary cases as well [30, 14]) shows that the steady-state current satisfies:

[TABLE]

If we use the notation $e_{i}:=\rho^{r}_{i}$, $\lambda_{0}:=\alpha$, and $\lambda_{n}:=\beta$ then this is just the steady-state equation of RFMEO given in (IV-D).