Optimal Translation Along a Circular mRNA

TL;DR

This paper introduces a spectral method for optimizing translation efficiency in circular mRNA models, enabling stable computation and convex optimization of ribosomal flow for genetic engineering.

Contribution

It derives a new spectral representation for the RFMR model, facilitating efficient and stable optimization of translation rates in circular mRNA.

Findings

Spectral representation simplifies optimal rate computation.

Optimal production rate is a strictly concave function.

Optimization becomes a convex problem with a unique solution.

Abstract

The ribosome flow model on a ring (RFMR) is a deterministic model for translation of a circularized mRNA. We derive a new spectral representation for the optimal steady-state production rate and the corresponding optimal steady-state ribosomal density in the RFMR. This representation has several important advantages. First, it provides a simple and numerically stable algorithm for determining the optimal values even in very long rings. Second, it enables efficient computation of the sensitivity of the optimal production rate to small changes in the transition rates along the mRNA. Third, it implies that the optimal steady-state production rate is a strictly concave function of the transition rates. Thus maximizing the optimal steady-state production rate with respect to the rates, under an affine constraint on the rates becomes a convex optimization problem that admits a unique solution, which can be determined numerically using highly efficient algorithms. This optimization problem is important, for example, when re-engineering heterologous genes in a host organism. We describe the implications of our results to this and other aspects of translation.

Figures

Equations

$${\dot{x}}_i(t)={\lambda }_{i-1}{x}_{i-1}(t)(1-x_i(t))-{\lambda }_i{x}_i(t)(1-x_{i+1}(t)),i=1,\dots ,n,$$

$$R:={\lambda }_n{e}_n.$$

$$B({\lambda }_0,\dots ,{\lambda }_n):=\left[\begin{array}{ccccccc}0& {\lambda }_0^{-1/2}& 0& 0& \dots & 0& 0\\ {\lambda }_0^{-1/2}& 0& {\lambda }_1^{-1/2}& 0& \dots & 0& 0\\ 0& {\lambda }_1^{-1/2}& 0& {\lambda }_2^{-1/2}& \dots & 0& 0\\ & & & ⋮& & & \\ 0& 0& 0& \dots & {\lambda }_{n-1}^{-1/2}& 0& {\lambda }_n^{-1/2}\\ 0& 0& 0& \dots & 0& {\lambda }_n^{-1/2}& 0\end{array}\right].$$

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

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

$$⋮$$

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

$$H(x(t)):=x_1(t)+\cdots +x_n(t),$$

$$L_s:=\{y\in C^n:y_1+\cdots +y_n=s\}.$$

$$\rho :=s/n$$

$${\lambda }_n{e}_n(1-e_1)$$

$$={\lambda }_2{e}_2(1-e_3),$$

$$⋮$$

$$={\lambda }_{n-1}{e}_{n-1}(1-e_n),$$

$$=R,$$

$$e_1+\cdots +e_n=s.$$

$$R(s,c{\lambda }_1,\dots ,c{\lambda }_n)=cR(s,{\lambda }_1,\dots ,{\lambda }_n),$$

$$\prod _{i=1}^n{e}_i^{\ast }=\prod _{i=1}^n(1-e_i^{\ast }).$$

$$e^{\ast }={\left[\begin{array}{ccc}0.6878& 0.3546& 0.4524\end{array}\right]}^{\prime },$$

$$A({\lambda }_1,\dots ,{\lambda }_n):=\left[\begin{array}{ccccccc}0& {\lambda }_1^{-1/2}& 0& 0& \dots & 0& {\lambda }_n^{-1/2}\\ {\lambda }_1^{-1/2}& 0& {\lambda }_2^{-1/2}& 0& \dots & 0& 0\\ 0& {\lambda }_2^{-1/2}& 0& {\lambda }_3^{-1/2}& \dots & 0& 0\\ & & & ⋮& & & \\ 0& 0& 0& \dots & {\lambda }_{n-2}^{-1/2}& 0& {\lambda }_{n-1}^{-1/2}\\ {\lambda }_n^{-1/2}& 0& 0& \dots & 0& {\lambda }_{n-1}^{-1/2}& 0\end{array}\right].$$

$$R^{\ast }$$

$$e_i^{\ast }$$

$$s^{\ast }$$

$$A=\left[\begin{array}{ccc}0& 1.1952& 0.6742\\ 1.1952& 0& 0.7906\\ 0.6742& 0.7906& 0\end{array}\right].$$

$$\zeta ={\left[\begin{array}{ccc}0.6024& 0.6219& 0.5004\end{array}\right]}^{\prime }.$$

$${\lambda }_1=\cdots ={\lambda }_n:={\lambda }_c,$$

$$\left[\begin{array}{ccccccc}0& {\lambda }_1^{-1/2}& 0& 0& \dots & 0& 0\\ {\lambda }_1^{-1/2}& 0& {\lambda }_2^{-1/2}& 0& \dots & 0& 0\\ 0& {\lambda }_2^{-1/2}& 0& {\lambda }_3^{-1/2}& \dots & 0& 0\\ & & & ⋮& & & \\ 0& 0& 0& \dots & {\lambda }_{n-2}^{-1/2}& 0& {\lambda }_{n-1}^{-1/2}\\ 0& 0& 0& \dots & 0& {\lambda }_{n-1}^{-1/2}& 0\end{array}\right].$$

$$e^{\ast }={\left[\begin{array}{ccccc}0.4260& 0.5831& 0.5939& 0.4019& 0.4950\end{array}\right]}^{\prime },R^{\ast }=\mathrm{0.1421.}$$

$$e^{\ast }={\left[\begin{array}{ccccc}0.9440& 0.7628& 0.6087& 0.2643& 0.0320\end{array}\right]}^{\prime },R^{\ast }=0.1791,$$

$$e^{\ast }={\left[\begin{array}{ccccc}0.9942& 0.7727& 0.6100& 0.2591& 0.0031\end{array}\right]}^{\prime },R^{\ast }=0.1808,$$

$$e^{\ast }={\left[\begin{array}{ccccc}0.9994& 0.7737& 0.6102& 0.2586& 0.0003\end{array}\right]}^{\prime },R^{\ast }=\mathrm{0.1810.}$$

$${\varphi }_i({\lambda }_1,\dots ,{\lambda }_n):=\frac{\partial }{\partial {\lambda }_i}{R}^{\ast }({\lambda }_1,\dots ,{\lambda }_n),i=1,\dots ,n,$$

$${\varphi }_i=\frac{2{\zeta }_i{\zeta }_{i+1}}{{\sigma }^3{\lambda }_i^{3/2}{\zeta }^{\prime }\zeta },i=1,\dots ,n.$$

$$\frac{{\varphi }_i}{{\varphi }_j}=\frac{{\zeta }_i{\zeta }_{i+1}}{{\zeta }_j{\zeta }_{j+1}}{\left(\frac{{\lambda }_j}{{\lambda }_i}\right)}^{3/2},i,j\in \{1,\dots ,n\},$$

$$\lambda ={\left[\begin{array}{ccccccc}1& 1.1& 0.55& 1.4& 1.3& 0.95& 0.6\end{array}\right]}^{\prime }.$$

Full text

I Introduction

Gene expression is the process by which the information encoded in a gene is used to synthesize a functional gene product. Two main stages of this process are transcription in which the information in the DNA of a specific gene is copied into a messenger RNA (mRNA) molecule and translation. The latter includes three phases: (1) initiation: complex macro-molecules called ribosomes bind to the mRNA;(2) elongation: the ribosomes unidirectionally decode each codon into the corresponding amino-acid that is delivered to the awaiting ribosome by transfer RNA (tRNA) molecules; and (3) termination: the ribosome detaches from the mRNA, the amino-acid sequence is released, folded and becomes a functional protein [3]. 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 occures in all living organisms, and under almost all conditions, to generate the macromolecular machinery for life. Developing a deeper understanding of translation has important implications in numerous scientific disciplines including medicine, evolutionary biology, biotechnology, and synthetic biology. Computational models of translation are essential in order to better understand 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., [12, 47, 46, 10, 41, 14, 36, 60]).

Computational models of translation describe the dynamical flow of ribosomes along the mRNA molecule, and include parameters that encode the various factors affecting the codon decoding rates and the binding of ribosomes. Some of these models provide a framework for both rigorous analysis and Monte Carlo simulations, thus promoting a better understanding of the way the parameters, and other factors, affect the dynamical and steady-state behavior of translation. Several computational models have been suggested based on different paradigms for example kinetics-based ordinary differential equations (see, e.g. [29]), Petri nets [7], and probabilistic Boolean networks [58]. For more details, see the survey papers [49, 60].

A standard mathematical model for ribosome flow is the totally asymmetric simple exclusion process (TASEP) [42, 59]. In this model, particles hop unidirectionally along an ordered lattice of $L$ sites. Each site can be either free or occupied by a particle, and a particle can only hop to a free site. This simple exclusion principle models particles that have “volume” and thus cannot overtake one other. The hops are stochastic, and the rate of hoping from site $i$ to site $i+1$ is denoted by $\gamma_{i}$. TASEP has two standard configurations. In TASEP with open boundary conditions the two sides of the lattice are connected to two particle reservoirs, and a particle can hop to [from] the first [last] site of the lattice at a rate $\alpha$ [$\beta$]. The average flow through the lattice converges to a steady-state value that depends on the parameters $\alpha,\gamma_{1},\dots,\gamma_{L-1},\beta$. Analysis of TASEP with open boundary conditions is non trivial, and closed-form results have been obtained mainly for the homogeneous TASEP (HTASEP), i.e. for the case where all the $\gamma_{i}$s are assumed to be equal.

In TASEP with periodic boundary conditions the chain is closed, and a particle that hops from the last site returns to the first one. Thus, here the lattice is a ring, and the total number of particles along the ring is conserved.

TASEP has become a fundamental model in non-equilibrium statistical mechanics, and has been applied to model numerous natural and artificial processes such as traffic flow, communication networks, and pedestrian dynamics [40]. In the context of translation, the lattice models the mRNA molecule, the particles are ribosomes, and simple exclusion means that a ribosome cannot overtake a ribosome in front of it.

The ribosome flow model (RFM) [39] is a continuous-time deterministic model for the unidirectional flow of “material” along an open chain of $n$ consecutive compartments (or sites). The RFM can be derived via a dynamic mean-field approximation of TASEP with open boundary conditions [40, section 4.9.7] [5, p. R345]. In a RFM with $n$ sites, the state variable $x_{i}(t)\in[0,1]$, $i=1,\dots,n$, describes the normalized amount of “material” (or density) at site $i$ at time $t$, where $x_{i}(t)=1$ [$x_{i}(t)=0$] indicates that site $i$ is completely full [completely empty] at time $t$. In the RFM, the two sides of the chain are connected to two particle reservoirs. A parameter $\lambda_{i}>0$, $i=0,\dots,n$, controls the transition rate from site $i$ to site $i+1$, where $\lambda_{0}$ [$\lambda_{n}$] is the initiation [exit] rate.

In the ribosome flow model on a ring (RFMR) [37] the particles exiting the last site reenter the first site. This is a dynamic mean-field approximation of TASEP with periodic boundary conditions. The RFMR admits a first integral, i.e. a quantity that is preserved along the dynamics, as the total ribosomal density is conserved. Both the RFM and RFMR are cooperative dynamical systems [43], but their dynamical properties are quite different [37].

Through simultaneous interactions with the cap-binding protein eIF4E and the poly(A)-binding protein PABP, the eukaryotic initiation factor eIF4G is able to bridge the two ends of the mRNA [52, 34]. This suggests that a large fraction of the ribosomes that complete translating the mRNA re-initiate. The RFMR is a good approximation of the translation dynamics in these circularized mRNAs. In addition, circular RNA forms (which includes covalent RNA interactions) appear in all domains of life [13, 11, 9, 8, 19, 1, 17, 2], and it was recently suggested that circular RNAs can be translated in eukaryotes [1, 17, 2].

It was shown in [37] that the RFMR admits a unique steady-state that depends on the total initial density and the transition rates, but not on the distribution of the total initial density among the sites. All trajectories emanating from initial conditions with the same total density converge to the unique steady-state. Ref. [56] considered the ribosomal density along a circular mRNA that maximizes the steady-state production rate using the RFMR. It was shown that given any arbitrary set of positive transition rates, there exists a unique optimal density (the same is true for TASEP with periodic boundary condition [28]). However, this unique optimum was not given explicitly, other than under certain special symmetry conditions on the rates, where the optimal density is one half of the maximal possible density.

The ribosomal density along the mRNA molecule plays a critical role in regulating gene expression, and specifically in determining protein production rates [30, 4]. For example, it was suggested in [4] that the cell tightly regulates ribosomal densities in order to maintain protein concentrations at different growth temperatures. At higher temperatures, the ribosomal density along the mRNA “improves” in order to increase protein production rates (as protein stability decreases with temperature).

The ribosomal density also affects different fundamental intracellular phenomena. Traffic james, abortions, and collisions may form if the ribosomal density is very high [44]. It may also contribute to co-translational misfolding of proteins, which then requires additional resources in order to degrade the degenerated proteins [15, 22, 57]. On the other hand, a very low ribosomal density may lead to high degradation rate of mRNA molecules [21, 16, 48, 35]. Thus, analyzing the ribosomal density that maximizes the production rate is critical in understanding how cells evolved to adapt and thrive in a changing environment.

Here we derive a new spectral representation (SR) for the optimal steady-state production rate and the corresponding steady-state ribosomal density in the RFMR. This SR has several important advantages. First, it provides a simple and numerically stable way to compute the optimal values even in very long rings. Second, it enables efficient computation of the sensitivity of the optimal steady-state production rate to small changes in the transition rates. This sensitivity analysis may find important applications in synthetic biology where a crucial problem is to determine the codons that are the most “important” in terms of their effect on the production rate. Third, the SR implies that the optimal steady-state production rate is a strictly concave function of the RFMR rates. Thus, the problem of maximizing the optimal steady-state production rate with respect to the rates becomes a convex optimization problem that admits a unique solution, which can be determined numerically using highly efficient algorithms.

The remainder of this paper is organized as follows. The next two sub-sections briefly review the RFM and the RFMR. Section II describes our main results and their biological implications. Section III concludes and suggests several directions for further research. To increase the readability of this paper, the proofs of all the results are placed in the Appendix.

I-A Ribosome Flow Model (RFM)

In a RFM with $n$ sites, the state variable $x_{i}(t)\in[0,1]$, $i=1,\dots,n$, denotes the density at site $i$ at time $t$, where $x_{i}(t)=1$ [$x_{i}(t)=0$] means that site $i$ is completely full [empty] at time $t$. The $n+1$ parameters $\lambda_{i}>0$, $i=0,\dots,n$, control the transition rate from site $i$ to site $i+1$. The RFM is a set of $n$ first-order nonlinear ordinary differential equations describing the change in the amount of “material” in each site:

[TABLE]

where $x_{0}(t):=1$, $x_{n+1}(t):=0$, and $\dot{x}_{i}$ is the change in the amount of material at site $i$ at time $t$, i.e. $\dot{x}_{i}(t):=\frac{d}{dt}x_{i}(t)$, $i=1,\dots,n$. Eq. (1) can be explained as a kind of a master equation: the change in density in site $i$ is the flow from site $i-1$ to site $i$ minus the flow from site $i$ to site $i+1$. The first flow, that is, the input rate to site $i$ is $\lambda_{i-1}x_{i-1}(t)(1-x_{i}(t))$. This rate is proportional to $x_{i-1}(t)$, i.e. it increases with the density at site $i-1$, and to $(1-x_{i}(t))$, i.e. it decreases as site $i$ becomes fuller. In particular, when site $i$ is completely full, i.e. when $x_{i}(t)=1$, there is no flow into this site. This is reminiscent of the simple exclusion principle: the flow of particles into a site decreases as that site becomes fuller. Note that the maximal possible flow from site $i-1$ to site $i$ is $\lambda_{i-1}$. Similarly, the output rate from site $i$, which is also the input rate to site $i+1$, is given by $\lambda_{i}x_{i}(t)(1-x_{i+1}(t))$. The output rate from the chain is $R(t):=\lambda_{n}x_{n}(t)$, that is, the flow out of the last site.

In the context of mRNA translation, the $n$-sites chain is the mRNA, $x_{i}(t)$ describes the ribosomal density at site $i$ at time $t$, and $R(t)$ describes the rate at which ribosomes leave the mRNA, which is also the rate at which the proteins are generated. Thus, $R(t)$ is the protein translation rate or production rate at time $t$.

Since every state-variable models the density of ribosomes in a site, normalized such that a value zero [one] corresponds to a completely empty [full] site, the state space of the RFM is the $n$-dimensional unit cube $C^{n}:=[0,1]^{n}$. Let $x(t,a)$ denote the solution of the RFM at time $t$ for the initial condition $x(0)=a$. It has been shown in [26] (see also [25]) that for every $a\in C^{n}$ this solution remains in $C^{n}$ for all $t\geq 0$, and that the RFM admits a globally asymptotically stable steady-state $e\in\operatorname{{\mathrm{i}nt}}(C^{n})$, i.e. $\lim_{t\to\infty}x(t,a)=e$ for all $a\in C^{n}$. The value $e$ depends on the rates $\lambda_{0},\dots,\lambda_{n}$, but not on the initial condition $x(0)=a$. This means that if we simulate the RFM starting from any initial density of ribosomes on the mRNA the dynamics will always converge to the same steady-state (i.e., to the same final ribosome density along the mRNA). In particular, the production rate $R(t)=\lambda_{n}x_{n}(t)$ always converges to the steady-state value:

[TABLE]

A spectral representation of this steady-state value has been derived in [32]. Given a RFM with dimension $n$ and rates $\lambda_{0},\dots,\lambda_{n}$, define a $(n+2)\times(n+2)$ Jacobi matrix

[TABLE]

Note that $B$ is componentwise non-negative and irreducible, so it admits a Perron root $\mu>0$. It has been shown in [32] that $\mu=R^{-1/2}$. This provides a way to compute the steady-state $R$ in the RFM without simulating the dynamical equations of the RFM.

For more on the analysis of the RFM using tools from systems and control theory and the biological implications of this analysis, see [54, 32, 33, 27, 53, 55]. Recently, a network of RFMs, interconnected via a pool of “free” ribosomes, has been used to model and analyze competition for ribosomes in the cell [38].

I-B Ribosome Flow Model on a Ring (RFMR)

If we consider the RFM under the additional assumption that all the ribosomes leaving site $n$ circulate back to site $1$ then we obtain the RFMR (see Fig. 2). Just like the RFM, the RFMR is described by $n$ nonlinear, first-order ordinary differential equations:

[TABLE]

The difference here with respect to the RFM is in the equations describing the change of material in sites $1$ and $n$. Specifically, the flow out of site $n$ is the flow into site $1$. This model assumes perfect recycling (be it covalent or non-covalent), and provides a good approximation when a large fraction of the ribosomes are recycled. Note that the RFMR can also be written succinctly as (1), but now with every index interpreted modulo $n$. In particular, $\lambda_{0}$ [$x_{0}$] is replaced by $\lambda_{n}$ [$x_{n}$].

Remark 1**.**

It is clear from the cyclic topology of the RFMR that if we cyclically shift all the rates $k$ times for some integer $k>1$ then the model does not change.

In the RFMR the total density of ribosomes along the ring at time $t$ is given by

[TABLE]

i.e. the sum of the density at each site. Let $s$ denote the total density of ribosomes along the chain at time $t=0$, i.e. $s:=H(x(0))$. Since ribosomes that exit site $n$ circulate back to site $1$, the total density is preserved for all time, that is, $H(x(t))\equiv s$ for all $t\geq 0$. The dynamics of the RFMR thus redistributes the particles between the sites, but without changing the total ribosome density. In the context of translation, this means that the total number of ribosomes on the (circular) mRNA is conserved. We say that $H(x(t))$ is a first integral of the RFMR.

For $s\in[0,n]$, the $s$ level set of $H$ is

[TABLE]

This is the set of all possible ribosome density configurations such that the total density is equal to $s$. For example, the vectors of densities $\begin{bmatrix}s&0&0&\dots&0\end{bmatrix}^{\prime}$ and $\begin{bmatrix}s/2&s/2&0&\dots&0\end{bmatrix}^{\prime}$ both belong to $L_{s}$.

Ref. [37] has shown that the RFMR is a strongly cooperative dynamical system, and that this implies that every level set $L_{s}$ contains a unique steady-state $e=e(s,\lambda_{1},\dots,\lambda_{n})\in\operatorname{{\mathrm{i}nt}}(C^{n})$, and that any trajectory of the RFMR emanating from any $x(0)\in L_{s}$ converges to this steady-state point. In particular, the production rate converges to a steady-state value $R=R(s,\lambda_{1},\dots,\lambda_{n})$.

Pick $s\in[0,n]$ and $a\in L_{s}$. Consider the RFMR with $x(0)=a$. Let

[TABLE]

denote the average ribosomal density in the RFMR. At steady-state, i.e. for $x=e$ the left-hand side of all the equations in (I-B) is zero, so

[TABLE]

and, since the total density is conserved,

[TABLE]

Note that it follows from (I-B) that for any $c>0$

[TABLE]

i.e. if we multiply all the rates by a factor $c>0$ then the steady-state production rate will also increase by the same factor $c$.

Given a set of transition rates, an interesting question is what ribosomal density maximizes the steady-state production rate in the RFMR? Indeed, a ribosomal density $s=0$ means zero production rate (as there are no ribosomes on the ring), and so does the completely full density $s=n$, as all the sites are completely full and the ribosomes cannot move forward. It was shown in [56] that for any arbitrary positive set of rates $\lambda_{1},\dots,\lambda_{n}$, there exists a unique density $s^{*}=s^{*}(\lambda_{1},\dots,\lambda_{n})$ (and thus a unique average density $\rho^{*}=s^{*}/n$) that maximizes the steady-state production rate. We denote the corresponding optimal steady-state production rate by $R^{*}=R(s^{*}(\lambda_{1},\dots,\lambda_{n}),\lambda_{1},\dots,\lambda_{n})$, and the corresponding optimal steady-state density by $e^{*}=e(s^{*}(\lambda_{1},\dots,\lambda_{n}),\lambda_{1},\dots,\lambda_{n})$. This means that in order to maximize the steady-state production rate (with respect to the total density), the mRNA must be initialized with a total density $s^{*}$ (the distribution of this total density along the mRNA at time zero is not important). Initializing with either more or less than $s^{*}$ (i.e with $\sum_{i=1}^{n}x_{i}(0)>s^{*}$ or with $\sum_{i=1}^{n}x_{i}(0)<s^{*}$ ) will decrease the steady-state production rate with respect to the one obtained when the circular mRNA is initialized with total density $\sum_{i=1}^{n}x_{i}(0)=s^{*}$.

The results in [56] also show that for the optimal value $s^{*}$, the steady-state density satisfies:

[TABLE]

This can be explained as follows. If the total density $s$ is too small then the product $\prod_{i=1}^{n}e_{i}$ is also small and thus $\prod_{i=1}^{n}e_{i}<\prod_{i=1}^{n}(1-e_{i})$. This case is not optimal i.e. it does not maximizes $R$, as there are not enough ribosomes on the ring. If the total density $s$ is too large then a similar argument yields $\prod_{i=1}^{n}e_{i}>\prod_{i=1}^{n}(1-e_{i})$. This case is also not optimal, as there are too many ribosomes on the ring and this leads to “traffic jams” that reduce the production rate. The optimal scenario lies between these two cases and is characterized by (7).

Example 1**.**

Fig. 3 depicts $R$ as a function of $s$ for a RFMR with dimension $n=3$ and rates $\lambda_{1}=0.7$, $\lambda_{2}=1.6$, and $\lambda_{3}=2.2$. It may be seen that there exists a unique value $s^{*}=1.4948$ (all numerical results in this paper are to four digit accuracy) that maximizes $R$. Simulating the RFMR with this initial density (e.g., by setting $x(0)=\begin{bmatrix}s^{*}/2&0&s^{*}/2\end{bmatrix}^{\prime}$) yields

[TABLE]

and $R^{*}=\lambda_{3}e^{*}_{3}(1-e^{*}_{1})=0.3107$. Note that $s^{*}$ is close (but not equal) to $3/2$, that is, one half of the maximal density. Note also that $\prod_{i=1}^{3}e^{*}_{i}=\prod_{i=1}^{3}(1-e^{*}_{i})=0.1103$.

Here, we present for the first time a spectral representation of the optimal steady-state production rate $R^{*}$ and the steady-state density $e^{*}$ in the non-homogeneous RFMR. We show that this representation has several advantages. First, it provides an efficient and numerically stable algorithm for evaluating $R^{*}$ and $e^{*}$ (and thus $s^{*}$) even for very large rings. This completely eliminates the need to simulate the RFMR dynamical equations for different values of $s$ in order to determine the optimal values. Furthermore, the spectral representation allows to analyze the sensitivity of $R^{*}$ to small changes in the rates. This sensitivity analysis could be crucial for example in synthetic biology applications, where an important problem is to determine positions along the transcript that affect the production rate the most (as this is not necessarily the positions of the slowest codons) [23]. Finally, we show that the spectral representation implies that $R^{*}$ is a strictly concave function of the rates. This means that the problem of maximizing $R^{*}$ with respect to the rates is a convex optimization problem, thus it admits a unique solution that can be efficiently determined numerically using algorithms that scale well with $n$.

It is important to note that in general the analysis results for the RFMR hold for any set of transition rates. This is in contrast to the analysis results for TASEP. Rigorous analysis of TASEP seems to be tractable only under the assumption that the internal hopping rates are all equal (i.e. the homogeneous case). In the context of translation, this models the very special case where all elongation rates are assumed to be equal.

The next section derives a spectral representation for $e^{*}$ and $R^{*}$, and describes its implications.

II Main Results

II-A Spectral Representation

Consider a RFMR with dimension $n>2$ and rates $\lambda_{1},\dots,\lambda_{n}>0$. Define an $n\times n$ matrix

[TABLE]

Note that this is a periodic Jacobi matrix (see, e.g. [51]).

We use the notation $\mathbb{R}^{n}_{++}:=\{v\in\mathbb{R}^{n}:v_{i}>0,i=1,\dots,n\}$, that is, the set of all $n$-dimensional vectors with positive entries. Since $A$ is symmetric, all its eigenvalues are real. Since $A$ is (componentwise) non-negative and irreducible, it admits a unique maximal eigenvalue $\sigma>0$ (called the Perron eigenvalue or Perron root), and a corresponding eigenvector $\zeta\in\mathbb{R}^{n}_{++}$ (the Perron eigenvector) [20].

Our first result provides a representation for the optimal steady-states in the RFMR using the spectral properties of the matrix $A$. In what follows, all indexes are interpreted modulo $n$. Recall that all the steady-state properties are invariant to any arbitrary cyclic shifts of the rates (see Remark 1), and that the proofs of all the results are placed in the Appendix.

Theorem 1**.**

Consider a RFMR with dimension $n$ and rates $\lambda_{1},\dots,\lambda_{n}$. Let $\sigma>0$ $[\zeta\in\mathbb{R}^{n}_{++}]$ denote the Perron eigenvalue [eigenvector] of $A$ in (8). Then the optimal values in the RFMR satisfy:

[TABLE]

Example 2**.**

Consider a RFMR with dimension $n=3$ and rates $\lambda_{1}=0.7$, $\lambda_{2}=1.6$, and $\lambda_{3}=2.2$. The corresponding matrix $A$ is:

[TABLE]

The maximal eigenvalue of $A$ is $\sigma=1.7940$, and the corresponding eigenvector is

[TABLE]

Now (1) yields $R^{*}=1.7940^{-2}=0.3107$, $e_{1}^{*}=0.6878$, $e_{2}^{*}=0.3546$, $e_{3}^{*}=0.4524$, and $s^{*}=1.4948$. This agrees of course with the results in Example 1.

Thm. 1 thus provides a spectral representation of the optimal values $R^{*}$, $e^{*}$, and $s^{*}$. This is important, since $e^{*}$ cannot be easily calculated based on the steady-state equations of the RFMR. Using simple and efficient algorithms to determine the eigenvalues and eigenvectors of a periodic Jacobi matrix, it is now possible to numerically compute $R^{*}$, $e^{*}_{i}$, $i=1,\dots,n$, and $s^{*}$ even for very large rings and without any simulations of the dynamical equations of the RFMR.

Thm. 1 has several more interesting implications. Given a RFMR with rates $\lambda_{1},\dots,\lambda_{n}$, define a vector $\bar{\lambda}\in\mathbb{R}^{n}_{++}$ by $\bar{\lambda}_{1}:=\lambda_{2}$, $\bar{\lambda}_{2}:=\lambda_{3},\dots,\bar{\lambda}_{n}=\lambda_{1}$. In other words, $\bar{\lambda}$ is a 1-step cyclic shift of $\lambda$. Let $P\in\mathbb{R}^{n\times n}$ be a matrix of zeros, except for the super-diagonal and the $(n,1)$ entry that are all equal to $1$. For example, for $n=4$, $P=\begin{bmatrix}0&1&0&0\\ 0&0&1&0\\ 0&0&0&1\\ 1&0&0&0\end{bmatrix}.$ Then $P$ is a permutation matrix so that $P^{\prime}=P^{-1}$, and $\bar{\lambda}=P\lambda$. It is straightforward to show that $A(\bar{\lambda})=PA(\lambda)P^{\prime}$, so $A(\lambda)$ and $A(\bar{\lambda})$ have the same spectral properties. Thus, Thm. 1 leads to the same steady-state results for both the original RFMR and its cyclic shift and this agrees with Remark 1.

In some special cases, the Perron eigenvalue and eigenvector of $A$ may be known explicitly and then one can immediately determine the optimal steady-state in the corresponding RFMR. The next example demonstrates this.

Example 3**.**

Consider a RFMR with homogeneous transition rates, i.e.

[TABLE]

where $\lambda_{c}$ denotes the common value of all the rates. Then it is straightforward to verify that $A(\lambda_{c},\cdots,\lambda_{c})$ admits a Perron eigenvalue $\sigma=2\lambda_{c}^{-1/2}$ and a corresponding eigenvector $\zeta=\begin{bmatrix}1&1&\cdots&1\end{bmatrix}^{\prime}$. Thm. 1 implies that $R^{*}=\lambda_{c}/4$ and $e_{i}^{*}=1/2$, $i=1,\dots,n$. This result has already been proven in [56, Prop. 3] using a different approach.

II-B Steady-State RFM as a Special Case of the Steady-State RFMR with Optimal Total Density

Comparing the spectral representations for the RFMR and the RFM yields the following result. Consider a RFMR with dimension $n$, fixed rates $\lambda_{1},\dots,\lambda_{n-1}$, and $\lambda_{n}\to\infty$. In this case, the matrix $A(\lambda)$ in (8) converges to the matrix:

[TABLE]

Comparing this with (3) and using Thm. 1 implies the following.

Corollary 1**.**

Let $e^{*}=\begin{bmatrix}e^{*}_{1}&\dots&e^{*}_{n}\end{bmatrix}^{\prime}$ denote the optimal steady-state of a RFMR with dimension $n$ and rates $\lambda_{1},\dots,\lambda_{n}$, where $\lambda_{n}\to\infty$. Let $\tilde{e}=\begin{bmatrix}\tilde{e}_{1}&\dots&\tilde{e}_{n-2}\end{bmatrix}^{\prime}$ denote the steady-state of a RFM with dimension $n-2$ and transition rates $\tilde{\lambda}_{0}=\lambda_{1},\tilde{\lambda}_{1}=\lambda_{2},\dots,\tilde{\lambda}_{n-2}=\lambda_{n-1}$. Then $\tilde{e}=\begin{bmatrix}e_{2}^{*}&e_{3}^{*}&\dots&e_{n-1}^{*}\end{bmatrix}^{\prime}$.

This implies that the steady-state of a RFM with arbitrary dimension $m$ and arbitrary rates $\tilde{\lambda}_{i}>0$ can be derived from the steady-state of an RFMR with dimension $n:=m+2$, rates $\lambda_{i}=\tilde{\lambda}_{i-1}$, $i=1,\dots,n-1$, $\lambda_{n}\to\infty$, that is initialized with the optimal total density $s^{*}$. In this respect, the RFM is a kind of “open-boundaries” RFMR that is initialized with the optimal total density.

This connection between the two models can be explained as follows. By (I-B), in an RFMR with $\lambda_{n}\to\infty$, the steady-state density at site $n$ will be zero, and at site $1$ it will be one. Indeed, the transition rate from site $n$ to site $1$ is infinite, so site $n$ will be completely emptied and site $1$ completely filled. This “disconnects” the ring at the link from site $n$ to site $1$. Furthermore, the completely full site $1$ serves as a “source” to site $2$ whereas the completely empty site $n$ serves as a “sink” to site $n-1$. The result is that sites $2,\dots,n-1$, of the RFMR become a RFM with dimension $n-2$. The next example demonstrates this.

Example 4**.**

Consider a RFMR with dimension $n=5$, and rates $\lambda_{1}=0.8,\lambda_{2}=0.6,\lambda_{3}=0.4,\lambda_{4}=0.7$, and $\lambda_{5}=0.5$. The optimal steady-state values are:

[TABLE]

For $\lambda_{5}=100$, the optimal steady-state values are:

[TABLE]

for $\lambda_{5}=10,000$, the optimal steady-state values are:

[TABLE]

and for $\lambda_{5}=1,000,000$, they are:

[TABLE]

It may be observed that as $\lambda_{5}$ increases, the optimal steady-state density at site $5$ [site $1$] decreases [increases] to zero [one]. On the other hand, for a RFM with dimension $n=3$ and rates $\tilde{\lambda}_{0}=0.8,\tilde{\lambda}_{1}=0.6,\tilde{\lambda}_{2}=0.4$, and $\tilde{\lambda}_{3}=0.7$, the steady-state values are: $\tilde{e}=\begin{bmatrix}0.7738&0.6102&0.2585\end{bmatrix}^{\prime}$, and $\tilde{R}=0.1810$ (compare to (12)).

II-C Sensitivity Analysis

We already know that given the transition rates $\lambda_{1},\dots,\lambda_{n}$, the RFMR admits a unique density $s^{*}(\lambda_{1},\dots,\lambda_{n})$ for which the steady-state production rate is maximized. Maximizing the steady-state production rate is a standard goal in biotechnology, and since codons may be replaced by their synonymous, an important question in the context of the RFMR is: how will a change in the rates affect the maximal production rate $R^{*}$? Note that the effect here is compound, as changing the rates also changes the optimal density that yields the maximal production rate.

In this section, we analyze

[TABLE]

i.e. the sensitivity of the optimal steady-state production rate $R^{*}$ with respect to $\lambda_{i}$.

A relatively large value of $\phi_{i}$ indicates that a small change in $\lambda_{i}$ will have a strong impact on the optimal steady-state production rate $R^{*}$. In other words, the sensitivities indicate which rates are the most “important” in terms of their effect on $R^{*}$. The results in Thm. 1 allow to compute the sensitivities using the spectral properties of the matrix $A$.

Proposition 1**.**

The sensitivities satisfy:

[TABLE]

Eq. (14) provides an efficient and numerically stable method to calculate the sensitivities for large-scale rings and arbitrary positive rates $\lambda_{i}$s using standard algorithms for computing the eigenvalues and eigenvectors of periodic Jacobi matrices. Note that (14) implies that all the sensitivities are positive.

Example 5**.**

Fig. 4 depicts $\ln(\phi_{i})$, computed using (14), as a function of $i$ for a RFMR with dimension $n=98$ and rates $\lambda_{1}=\lambda_{50}=0.3$ and $\lambda_{i}=1$ for all other $i$. Here the maximal sensitivity is $\phi_{1}=\phi_{50}$, and the sensitivities decrease as we move away from sites $1$ and $50$. This makes sense as the corresponding rates are the bottleneck rates in this example.

Eq. (14) implies that

[TABLE]

that is, the ratio between any two sensitivities is determined by the corresponding Perron eigenvalue components and the corresponding rates. One may expect that the highest sensitivity will correspond to the minimal rate, but (15) shows that this is not necessarily so. The next example demonstrates this.

Example 6**.**

Consider a RFMR with dimension $n=7$ and rates:

[TABLE]

In this case, $R^{*}=0.2213$. Using (14) yields the sensitivities:

[TABLE]

Note that although the minimum rate is $\lambda_{3}$, the maximal sensitivity is $\phi_{7}$. This implies that increasing $\lambda_{7}$ by some small value $\varepsilon>0$ will increase $R^{*}$ more than the increase due to increasing any other rate by $\varepsilon$. For example, increasing $\lambda_{3}$ by $0.05$ (and leaving all other rates unchanged) yields $R^{*}=0.2248$, while increasing $\lambda_{7}$ by $0.05$ instead (and leaving all other rates unchanged) yields $R^{*}=0.2251$.

The spectral approach can also be used to derive theoretical results on the sensitivities. The next three results demonstrate this.

Proposition 2**.**

*The sensitivities satisfy $0<\phi_{i}\leq 1$ for all $i=1,\dots,n$. *

This implies that an increase [decrease] in any of the rates by $\varepsilon$ increases [decreases] the optimal steady-state production rate by no more than $\varepsilon$.

Proposition 3**.**

Consider a RFMR with dimension $n$ and homogeneous rates (10). Then

[TABLE]

This means that in the homogeneous case, all the sensitivities are equal. This is of course expected, as the circular topology of the sites implies that all the rates have the same effect on $R^{*}$. Furthermore, the sensitivities decrease with $n$, i.e. in a longer ring each rate has a smaller effect on $R^{*}$.

Assume now that the RFMR rates satisfy

[TABLE]

i.e. the rates are symmetric. Note that since all indexes are interpreted modulo $n$, it is enough that (16) holds for some cyclic permutation of the rates. For example, for $n=3$ the rates are symmetric if at least two of the rates $\lambda_{1},\lambda_{2},\lambda_{3}$ are equal.

Proposition 4**.**

Consider a RFMR with dimension $n$ and symmetric rates (16). Then

[TABLE]

Example 7**.**

Consider a RFMR with dimension $n=6$ and rates $\lambda_{1}=\lambda_{5}=1$, $\lambda_{2}=\lambda_{4}=1.2$, $\lambda_{3}=0.8$ and $\lambda_{6}=1.5$. Note that these rates satisfy (16). The sensitivities are:

[TABLE]

and it may be observed that $\phi_{i}=\phi_{6-i}$, $i=1,\dots,5$.

II-D Optimizing the production rate

Any set of rates $\lambda=(\lambda_{1},\dots,\lambda_{n})$ induces an optimal density $s^{*}$ and the RFMR initialized with this total density yields a maximal production rate $R^{*}$ (with respect to all other initial densities). This yields a mapping $\lambda\to R^{*}(\lambda)$. Now suppose that we have some set, denoted by $\Omega$, of $n$-dimensional vectors with positive entries. Every vector from $\Omega$ can be used as a set of rates $\lambda$ for the RFMR, and thus yields a value $R=R^{*}(\lambda)$. A natural question is: determine a vector $\eta\in\Omega$ that yields the maximal value, that is,

[TABLE]

In the context of translation, this means that a circular mRNA with rates $\eta$, initialized with total density $s^{*}(\eta)$, will yield a steady-state production rate that is higher than that obtained for all the other options for the rate vector in $\Omega$ (regardless of the initial total density in these other circular mRNAs).

The next result is essential for efficiently analyzing the maximization of $R^{*}$ with respect to (w.r.t.) its rates.

Proposition 5**.**

Consider a RFMR with dimension $n$. The mapping $\lambda:=(\lambda_{1},\dots,\lambda_{n})\mapsto R^{*}(\lambda)$ is strictly concave on $\mathbb{R}^{n}_{++}$.

For example, for $n=2$ it is straightforward to show that

[TABLE]

Fig. 5 depicts $R^{*}(\lambda_{1},\lambda_{2})$ as a function of its parameters. It may be observed that this is a strictly concave function on $\mathbb{R}^{2}_{++}$.

The sensitivity analysis of $R^{*}$, and its strict concavity w.r.t the rates, have important implications to the problem of optimizing the steady-state production rate in the RFMR w.r.t the rates $\lambda$. We now explain this using a specific optimization problem. First note that to make the problem meaningful every rate must be bounded above. Otherwise, the optimal solution will be to take this rate to infinity. We thus consider the following constrained optimization problem.

Problem 1**.**

Consider a RFMR with dimension $n$. Given the parameters $w_{1},\dots,w_{n},b>0$, maximize $R^{*}=R^{*}(\lambda_{1},\dots,\lambda_{n})$ with respect to the parameters $\lambda_{1},\dots,\lambda_{n}$, subject to the constraints

[TABLE]

In other words, the problem is to maximize $R^{*}$ w.r.t. the rates, under the constraints that the rates are positive and their weighted sum is bounded by $b$. The weights $w_{i}$s can be used to provide different weighting to the different rates, and $b$ represents the “total biocellular budget”. By Prop. 2, the optimal solution always satisfies the constraint in (17) with equality. Note that a similar optimization problem was defined and analyzed in the context of the RFM in [32].

In the context of mRNA translation, each $\lambda_{i}$ depends on the availability of translation resources that affect codon decoding times, such as tRNA molecules, amino acids, elongation factors, and Aminoacyl tRNA synthetases. These resources are limited as generating them consumes significant amounts of cellular energy. They are also correlated. For example, a large $\lambda_{i}$ may imply large consumption of certain tRNA molecules by site $i$, depleting the availability of tRNA molecules to the other sites. Thus, the first (affine) constraint in (17) describes the limited and shared translation resources, whereas $b$ describes the total available biocellular budget.

By Prop. 5, the objective function in Problem 1 is strictly concave, and since the constraints are affine, Problem 1 is a convex optimization problem [6]. Thus, it admits a unique solution. We denote the optimal solution of Problem 1 by $\lambda^{co}:=(\lambda_{1}^{co},\dots,\lambda_{n}^{co})$, and the corresponding maximal (now in the sense of total density and transition rates) steady-state production rate by $R^{co}$ (where $co$ denotes constrained optimization). This means that for a RFMR with dimension $n$, $R^{co}$ is the maximal steady-state production rate over all the rates satisfying the constraints (17) and all possible total initial densities.

The convexity also implies that the solution can be determined efficiently using numerical algorithms that scale well with $n$. To demonstrate this, we wrote a simple and unoptimized MATLAB program (that is guaranteed to converge because of the convexity) for solving this optimization problem and ran it on a MAC laptop with a $2.6$ GHz Intel core i7 processor. As an example, for $n=100$ and the (arbitrarily chosen) weights $w_{i}=1+0.4\sin(2\pi i/100)$, $i=1,\dots,100$, and $b=1$, the optimal solution was found after $11.7$ seconds.

The affine constraint in (17) includes a possibly different weight for each of the rates. For example, if $w_{2}$ is much larger than the other weights then this means that any small increase in $\lambda_{2}$ will greatly increase the total weighted sum, thus typically forcing the optimal value $\lambda_{2}^{co}$ to be small. In the special case where all the $w_{i}$s are equal the formulation gives equal preference to all the rates, so if the corresponding optimal solution satisfies $\lambda_{j}^{co}>\lambda_{i}^{co}$, for some $i,j$, then this implies that, in the context of maximizing $R^{*}$, $\lambda_{j}$ is “more important” than $\lambda_{i}$. We refer to this case as the homogeneous constraint case and assume, without loss of generality, that $w_{i}=1$ for all $i$. Note that by (6) we can always assume, without loss of generality, that $b=1$.

Proposition 6**.**

Consider Problem 1 with $w_{1}=\cdots=w_{n}=b=1$, i.e. the affine constraint is

[TABLE]

Then the optimal solution is $\lambda_{i}^{co}=1/n$ for all $i$. The RFMR with these rates satisfies $s^{*}=n/2$, $e_{i}^{co}=1/2$ for all $i$, and $R^{co}=1/(4n)$.

Remark 2**.**

In view of the Kuhn–Tucker theorem [6], the necessary and sufficient condition for optimality of $\lambda$ in Problem 1 with homogeneous weights is that the sensitivity $\phi_{i}=\frac{\partial R^{*}}{\partial\lambda_{i}}(\lambda^{*})$ does not depend on the index $i$.

III Discussion

We considered a deterministic model for translation along a circular mRNA. The behavior of this model depends on the transition rates between the sites and on the value $s:=\sum_{i=1}^{n}x_{i}(0)$, that is, the initial total density along the ring. The total density is conserved, so $\sum_{i=1}^{n}x_{i}(t)=s$ for all $t\geq 0$.

We derived a spectral representation for the steady-state density and production rate for the case where the initial density is $s^{*}$, i.e. the density yielding a maximal steady-state production rate. In fact, the proof of Thm. 1 (see the Appendix) shows that we can interpret the optimal density RFMR as a dynamical system that “finds” the Perron eigenvalue and eigenvector of a certain periodic Jacobi matrix.

The spectral representation for the RFMR provides a powerful framework for analyzing the RFMR when initialized with the optimal total density $s^{*}$. In addition to providing an efficient and numerically stable manner for computing the optimal steady-state production rate and steady-state density, it allows to efficiently compute the sensitivity of the optimal steady-state production rate to perturbations in the rates. This is important as conditions in the cell are inherently stochastic, and thus sensitivity analysis must accompany the steady-state description.

Furthermore, using the spectral representation, it was shown that the steady-state production rate with optimal density is a strictly concave function of the RFMR rates. The translation machinery in the cell is affected by different kinds of mutations (e.g. synonymous codon mutations, duplication of a tRNA gene, etc.). The strict concavity result thus suggest that the selection of mutations that increase fitness indeed converges towards the unique optimal parameter values (by a simple “hill-climbing” evolution process). The strict concavity implies that given an affine (and more generally convex) constraint on the rates, that represents limited and shared translation resources, the unique optimal set of rates can be determined efficiently even for (circular) mRNAs with a large number of codons.

Obtaining an optimal production rate is an important problem in synthetic biology and biotechnology. Examples include optimal synonymous codon mutations of an endogenous gene, and optimal translation efficiency and protein levels of heterologous genes in a new host [31, 45, 18, 21]. These genes compete with endogenous genes for the available translation resources, as consuming too much resources by the heterologous gene may kill the host [31, 45]. Thus, any realistic optimization of the protein production rate should not consume too many resources, as otherwise the fitness of the host may be significantly reduced. These considerations seems to fit well with the affine-constrained optimization problem presented and analyzed here.

We also showed that the spectral representation of the RFM follows as a special case of the representation given here for the RFMR. However, it seems that a better understanding of the link between the RFM and the RFMR requires further study. Our results suggest several other interesting directions for future research. One such direction is finding special cases, besides the one described in Example 3, where the Perron eigenvalue and eigenvector of $A(\lambda_{1},\dots,\lambda_{n})$ are explicitly known. Another possible direction is the analysis of the dual of the optimization problem defined by Problem 1. Specifically, does the dual problem has any interesting biological interpretation in the context of mRNA translation, and does its analysis provides more insight into optimizing translation?

Finally, TASEP with periodic boundary conditions has been used to model many transport phenomena including traffic flow and pedestrian dynamics [40, 50]. We believe that the spectral representation of the RFMR with optimal density may be useful also for analyzing these transport applications.

Acknowledgments

The research of YZ is partially supported by the Edmond J. Safra Center for Bioinformatics at Tel Aviv University. The research of AO is partially supported by the Russian Foundation for Basic Research, grant 17-08-00742. The research of MM is partially supported by research grants from the Israeli Ministry of Science, Technology & Space, the US-Israel Binational Science Foundation, and the Israeli Science Foundation.

Author Contributions Statement

YZ, AO, and MM performed the research and wrote the paper.

Data Availability Statement

All the relevant data is included in the manuscript.

Competing Financial Interests Statement

The authors declare no competing financial interests.

Appendix - Proofs

• Proof of Thm. 1.

Pick $n>2$ and parameters $c_{1},\dots,c_{n-1}>0$, and $c_{n}\geq 0$. Consider the $n\times n$ periodic Jacobi matrix:

[TABLE]

Note that $J$ is irreducible and (componentwise) non-negative. Let $\sigma>0$ denote that Perron eigenvalue of $J$ and let $\zeta\in\mathbb{R}^{n}_{++}$ denote the corresponding eigenvector. The equation $J\zeta=\sigma\zeta$ yields

[TABLE]

Define

[TABLE]

Note that since the indexes are interpreted modulo $n$, Eq. (20) implies in particular that

[TABLE]

Then (Proof of Thm. 1.) yields:

[TABLE]

Also, it follows from (20) that $\prod_{i=1}^{n}d_{i}=\sigma^{-n}\prod_{i=1}^{n}c_{i}$, and from (Proof of Thm. 1.) that $\prod_{i=1}^{n}(1-d_{i})=\sigma^{-2n}\frac{\prod_{i=1}^{n}c_{i}^{2}}{\prod_{i=1}^{n}d_{i}}$, and combining these two equations yields

[TABLE]

Note that all the derivations above hold for any real eigenvalue of $J$ and its corresponding eigenvector (assuming all its entries are non zero so that (20) is well-defined), but since the Perron eigenvector is the only eigenvector in the first orthant [20], all the $d_{i}$s are positive only for the Perron eigenvalue and eigenvector.

Now consider a RFMR with dimension $n$ and rates $\lambda_{i}:=c_{i}^{-2}$, $i=1,\dots,n$, that is:

[TABLE]

We already know that this system converges to a steady-state $e\in C^{n}$, that is,

[TABLE]

Comparing this with (Proof of Thm. 1.) shows that $e_{i}=d_{i}$ for all $i$, and that the steady-state production rate is $R=\sigma^{-2}$. Furthermore, (23) implies that $\prod_{i=1}^{n}e_{i}=\prod_{i=1}^{n}(1-e_{i})$, so we conclude that the steady-state satisfies condition (7) that describes the unique optimal steady-state (i.e. the steady-state production rate that corresponds to the unique optimal total density $s^{*}$). This proves the first two equations in (1). Finally, since the total density is conserved, it is equal to $\sum_{i=1}^{n}e_{i}$. This completes the proof of Thm. 1. ∎

• Proof of Proposition 1.

By Thm. 1,

[TABLE]

By known results from linear algebra (see, e.g., [24]), the sensitivity of the Perron root of $A$ with respect to a change in $\lambda_{i}$ is

[TABLE]

Only the entries $a_{i,i+1}=a_{i+1,i}=\lambda_{i}^{-1/2}$ depend on $\lambda_{i}$, so

[TABLE]

and combining this with (25) proves (14). ∎

• Proof of Prop. 2.

Since $\sigma>0$ and $\zeta\in\mathbb{R}^{n}_{++}$, $\phi_{i}>0$ for all $i$. To prove the upper bound, perturb $\lambda_{i}$ to $\bar{\lambda}_{i}:=\lambda_{i}+\varepsilon$, with $\varepsilon>0$ sufficiently small. This yields a perturbed matrix $\bar{A}$ that is identical to $A$ except for entries $(i,i+1)$ and $(i+1,i)$ that are

[TABLE]

where $o(\varepsilon)$ denotes a function $f(\varepsilon)$ satisfying $\lim_{\varepsilon\to 0}\frac{f(\varepsilon)}{\varepsilon}=0$. This means that $\bar{A}=A+P$, where $P\in\mathbb{R}^{n\times n}$ is a matrix with zero entries except for entries $(i,i+1)$ and $(i+1,i)$ that are equal to $-\frac{\varepsilon\lambda_{i}^{-3/2}}{2}+o(\varepsilon)$. By Weyl’s inequality [20], $\rho(\bar{A})\geq\rho(A)-\frac{\varepsilon\lambda_{i}^{-3/2}}{2}+o(\varepsilon)$, where $\rho(Q)$ denotes the maximal eigenvalue of a symmetric matrix $Q$. This means that $\frac{\partial\rho(A)}{\partial\lambda_{i}}\geq-\frac{\lambda_{i}^{-3/2}}{2}+\frac{o(\varepsilon)}{\varepsilon}$, thus $\phi_{i}\leq(R^{*}/\lambda_{i})^{3/2}$. Since $R^{*}\leq\lambda_{i}$, it follows that $\phi_{i}\leq 1$ for all $i$. ∎

• Proof of Prop. 3.

Consider a RFMR with homogeneous rates (10). Then by Example 3, $\zeta_{i}=1$, $i=1,\dots,n$, and $\sigma=2\lambda_{c}^{-1/2}$, and plugging these in (14) completes the proof.

∎

• Proof of Prop. 4.

We require the following result.

Proposition 7**.**

Consider the RFMR with dimension $n$ and symmetric rates. Then $\zeta_{i}=\zeta_{n+1-i}$, $i=1,\dots,n$.

• Proof of Prop. 7.

Consider first the case $n$ even. Let $Q\in\mathbb{R}^{(n/2)\times(n/2)}$ be a reversal matrix, i.e. a matrix of zeros except for the counter-diagonal (i.e. entries $(i,\frac{n}{2}-i+1)$, $i=1,\dots,n/2$) that is all ones. For example, for $n=4$,

[TABLE]

Note that given any arbitrary vector $v=\begin{bmatrix}v_{1}&v_{2}&\cdots&v_{n/2}\end{bmatrix}^{\prime}\in\mathbb{R}^{n/2}$, $Qv=\begin{bmatrix}v_{n/2}&v_{(n/2)-1}&\cdots&v_{1}\end{bmatrix}^{\prime}$.

Since the rates satisfy (16), the $n\times n$ matrix $A$ has the form

[TABLE]

where $A_{1}\in\mathbb{R}^{(n/2)\times(n/2)}_{+}$ is a matrix of zeros except for the super-diagonal and the sub-diagonal, which are both equal to $(\lambda_{1}^{-1/2},\dots,\lambda_{(n/2)-1}^{-1/2})$, and $A_{2}\in\mathbb{R}^{(n/2)\times(n/2)}_{+}$ is a matrix of zeros except for entry $(1,n/2)$ that is $\lambda_{n}^{-1/2}$, and entry $(n/2,1)$ that is $\lambda_{n/2}^{-1/2}$. Decompose the Perron eigenvector $\zeta$ of $A$ as $\zeta^{1}:=\begin{bmatrix}\zeta_{1}&\dots&\zeta_{n/2}\end{bmatrix}^{\prime}$ and $\zeta^{2}:=\begin{bmatrix}\zeta_{(n/2)+1}&\dots&\zeta_{n}\end{bmatrix}^{\prime}$.

Let $\rho(W)$ denote the spectral radius of a matrix $W$. Since $A_{1}$ is a principal submatrix of the (componentwise) nonnegative matrix $A$, $\rho(A_{1})\leq\rho(A)$ (see [20, Ch. 8]). Assume for the moment that $\rho(A_{1})=\rho(A)$. Then using the fact that $QQ=I$, that is $Q=Q^{-1}$, we conclude that $\rho\left(\begin{bmatrix}A_{1}&0\\ 0&QA_{1}Q\end{bmatrix}\right)=\rho(A)$. This means that the matrices $\begin{bmatrix}A_{1}&0\\ 0&QA_{1}Q\end{bmatrix}$ and $\begin{bmatrix}A_{1}&A_{2}\\ QA_{2}Q&QA_{1}Q\end{bmatrix}$ have the same Perron root, but this contradicts Prop. 2. We conclude that

[TABLE]

The equation $A\zeta=\sigma\zeta$ yields

[TABLE]

Multiplying both sides of the second equation by $Q$, noting that $QQ=I$, and rearranging yield

[TABLE]

Subtracting the second equation from the first and using again the fact that $QQ=I$ yields

[TABLE]

Combining this with (26) and the fact that $A_{2}Q$ is (componentwise) nonnegative implies that $\zeta^{1}=Q\zeta^{2}$, i.e. $\zeta_{i}=\zeta_{n+1-i}$, $i=1,\dots,n$. This completes the proof for the case $n$ even. The proof when $n$ is odd is very similar and therefore omitted. ∎

Now the proof of Prop. 4 follows from combining (14), Thm. 1, and Prop. 7. ∎

• Proof of Prop. 5.

Indeed, the map $\lambda=(\lambda_{1},\dots,\lambda_{n})\mapsto A(\lambda)$, where $A(\lambda)$ is given in (8), from $\mathbb{R}^{n}_{++}$ is convex, meaning that the matrix inequality

[TABLE]

holds elementwise for any arbitrary $\lambda^{\prime},\lambda^{\prime\prime}\in\mathbb{R}^{n}_{++}$. This immediately follows from the convexity of the real function $\lambda\mapsto\lambda^{-1/2}$. The Perron–Frobenius theorem implies the corresponding inequality for the Perron eigenvalue [20]

[TABLE]

where the inequality (30) is strict if $\lambda^{\prime}\neq\lambda^{\prime\prime}$. Thus $\sigma(\lambda)$ is a strictly convex function. In view of the basic identity $R^{*}=\sigma^{-2}$ in (1), it follows that $R^{*}(\lambda)$ is a strictly concave function. ∎

• Proof of Prop. 6.

We know that Problem 1 admits a unique optimal solution $\tilde{\lambda}$. Consider the cyclic shift $\bar{\lambda}_{i}=\tilde{\lambda}_{i+1}$, $i=1,\dots,n$, where the indices are taken modulo $n$. Note that $\sum_{i=1}^{n}\bar{\lambda}_{i}=\sum_{i=1}^{n}\tilde{\lambda}_{i}=1$, so $\bar{\lambda}$ also satisfies the constraint (18). The matrices $A(\bar{\lambda})$ and $A(\tilde{\lambda})$ have the same spectrum. Since the optimal solution is unique, $\overline{\lambda}=\tilde{\lambda}$. We conclude that the optimal transition rates $\tilde{\lambda}_{i}$ are all equal, and thus $\lambda_{i}^{co}:=1/n$, $i=1,\dots,n$. By Example 3, $R^{co}=1/(4n)$, and $e^{co}_{i}=1/2$, $i=1,\dots,n$. ∎