Translation in the cell under fierce competition for shared resources: a mathematical model

During translation, mRNAs ‘compete’ for shared resources. Under stress conditions, during viral infection and also in high-throughput heterologous gene expression, these resources may become scarce, e.g. the pool of free ribosomes is starved, and then the competition may have a dramatic effect on the global dynamics of translation in the cell. We model this scenario using a network that includes m ribosome flow models (RFMs) interconnected via a pool of free ribosomes. Each RFM models ribosome flow along an mRNA molecule, and the pool models the shared resource. We assume that the number of mRNAs is large, so many ribosomes are attached to the mRNAs, and the pool is starved. Our analysis shows that adding an mRNA has an intricate effect on the total protein production. The new mRNA produces new proteins, but the other mRNAs produce less proteins, as the pool that feeds these mRNAs now has a smaller abundance of ribosomes. As the number of mRNAs increases, the marginal utility of adding another mRNA diminishes, and the total protein production rate saturates to a limiting value. We demonstrate our approach using an example of insulin protein production in a cell-free system.

TT, 0000-0003-4194-7068; MM, 0000-0001-8319-8996 During translation, mRNAs 'compete' for shared resources. Under stress conditions, during viral infection and also in high-throughput heterologous gene expression, these resources may become scarce, e.g. the pool of free ribosomes is starved, and then the competition may have a dramatic effect on the global dynamics of translation in the cell. We model this scenario using a network that includes m ribosome flow models (RFMs) interconnected via a pool of free ribosomes. Each RFM models ribosome flow along an mRNA molecule, and the pool models the shared resource. We assume that the number of mRNAs is large, so many ribosomes are attached to the mRNAs, and the pool is starved. Our analysis shows that adding an mRNA has an intricate effect on the total protein production. The new mRNA produces new proteins, but the other mRNAs produce less proteins, as the pool that feeds these mRNAs now has a smaller abundance of ribosomes. As the number of mRNAs increases, the marginal utility of adding another mRNA diminishes, and the total protein production rate saturates to a limiting value. We demonstrate our approach using an example of insulin protein production in a cell-free system.

Introduction
mRNA translation is a fundamental process in gene expression, i.e. the transformation of genetic information into a functional protein [1]. During translation, the ribosomes (complex macro-molecules) scan the mRNA in a sequential manner, 'read' it codon by codon, and build a corresponding chain of amino acids. When the ribosome reaches a stop codon, it detaches from the mRNA and releases the chain that, after additional processing, becomes a functional protein.
Many ribosomes may attach to the same mRNA and translate it in parallel. This pipelining increases the protein production rate. The speed of the ribosome movement along the mRNA is determined by the mRNA sequence and structure, pauses due to collisions with other ribosomes, and by various translation factors, e.g. the abundance of cognate tRNA molecules in the vicinity of the mRNA. Understanding the dynamics of translation is of considerable importance, as it plays a major role in determining the protein production rate [2,3]. Furthermore, the dynamics of mRNA translation and specifically the evolution of ribosomal 'traffic jams' [4,5] along an mRNA have been implicated to various diseases [6,7].
Computational models for ribosome flow along an mRNA and the resulting protein production rate can be used to integrate and explain the growing number of experimental findings (e.g. via methods like ribosome profiling [8] that can be performed for single cells [9], and methods for imaging the translation of a single mRNA molecule [10]). Such models can also predict the effect of various manipulations or regulation of the genetic machinery on the protein production rate. Such manipulations are common both in biotechnology and also during viral infection, where the virus 'hijacks' and potentially shuts down parts of the host cell translation machinery. For example, SARS-CoV-2 uses a multipronged strategy to impede host protein synthesis and affect translation in a global manner [11]. In addition, heterologous genes tend to overload the translation machinery and also have a global effect on translation [12,13].
A popular model for ribosome flow, and many other natural and artificial systems and processes, is the totally asymmetric simple exclusion process (TASEP) (see e.g. [14,15]). This is a stochastic discrete-time process that includes a onedimensional chain of sites and particles that hop stochastically along the chain in a unidirectional manner. Simple exclusion refers to the fact that two particles cannot be in the same site at the same time, i.e. a particle can only hop to an empty site. In the context of translation, the particles are ribosomes and the sites are consecutive (groups of) codons along the mRNA [16][17][18]. Simple exclusion corresponds to the fact that a ribosome cannot 'overtake' a ribosome in front of it, as the information on the mRNA must be decoded in a sequential manner. TASEP has become a phenomenological model in statistical mechanics, yet rigorous analysis of this model is difficult, except for some very special cases, e.g. when all the internal hopping rates are assumed to be equal [14].
The ribosome flow model (RFM) is the dynamic mean-field approximation of TASEP. This yields a continuous-time, deterministic, nonlinear model for the flow of ribosomes along the mRNA [19]. The RFM is highly amenable to analysis using tools from systems and control theory, and it has been used to model and analyse many important aspects of translation including: entrainment of the protein production rate to periodic initiation and elongation rates with a common period [20], sensitivity analysis of the steady-state production rate [21], optimizing the protein production rate subject to convex constraints on the rates [22], the effect of ribosome recycling [23,24], determining the ribosome density that maximizes protein production [25], stochastic variability in translation [26], maximizing the average protein production [27] and more [28].
The cell is a factory for producing proteins that includes a large number of mRNA molecules and ribosomes. For example, a Saccharomyces cerevisiae cell includes about 60 000 mRNA molecules and 200 000 ribosomes. About 85% of the ribosomes are associated with mRNAs [29][30][31], and the rest form the 'pool of free ribosomes'. The mRNAs thus indirectly 'compete' for the available ribosomes. This generates a network with intricate indirect coupling between the mRNAs. For example, ribosomal 'traffic jams' on mRNAs may deplete the pool leading to lower initiation rates in the other mRNAs. When the shared resources are abundant the coupling between mRNAs due to this competition is weak, and the network can potentially be analysed using models of translation along a single, isolated RNA. However, when the resources are scarce, e.g. when the pool of free ribosomes is starved, the competition may have a strong effect on the global dynamics of translation in the cell.
The latter scenario may be relevant both in synthetic systems, where the goal is to optimize the production rate, and under various physiological conditions. For example, under stress conditions or during a high-yield viral infection, where the viral mRNAs 'hijack' the translation machinery, and consume many of the shared resources to produce viral proteins. It may also be relevant in heterologous gene expression (e.g. when the heterologous gene is overexpressed and consumes most of the ribosomes in the cell), and in cell-free systems where the number of mRNA molecules may be relatively large in comparison with the number of ribosomes [32].
We study this scenario using a mathematical model that includes a network of m RFMs interconnected via a pool of free ribosomes. The RFM belongs to the class of compartmental models that have been used in various domains of systems biology including genetics, physiology and pharmacology [33]. Each RFM describes the dynamics of translation along one mRNA molecule, and the interconnection via the pool encapsulates the competition for shared resources. This model was first suggested in [34]. Mathematically, it is a cooperative dynamical system [35,36] that admits a first integral: the total density of ribosomes in the network is conserved. Such systems have a well-ordered asymptotic behaviour [37,38]. It was shown in [34] that any solution of the network converges to a steady state, where the flow of ribosomes into and out of any site along any mRNA is equal. Also, the flows in and out of the pool are equal. Sensitivity analysis of this steady state [34] with respect to modifying one of the translation rates in a specific mRNA demonstrated that the steady-state production rates in all the other mRNAs either all increase or all decrease.
The model based on networks of RFMs has already been validated experimentally, and applied to predict the density of ribosomes along different mRNAs, the protein levels of different genes, and even for re-engineering ribosomal traffic jams (see e.g. [19,31]).
A generalization of this network, that includes the possibility of ribosome drop-off and attachment at each site along the RFM, was recently suggested in [39]. Simulations of this model showed that ribosome drop-off from an isolated mRNA always decreases the protein production rate, yet in the network ribosome drop-off from a jammed mRNA may increase the total production rate of all the mRNAs, as the drop-off frees ribosomes that enter the pool, and this improves the production rate in the other mRNAs. This illustrates how the network perspective provides new biological insights.
In this paper, we analyse this network of m RFMs interconnected via a pool from a new, structural perspective, that is, we study how adding new RFMs to the network affects the dynamics. The main contributions of this paper include the following: (1) We prove that adding an RFM to the network always decreases the steady-state pool density. This makes sense, as every new RFM 'consumes' ribosomes. (2) We show that adding an RFM has an intricate effect on the network steady state: on the one hand, the additional RFM produces new proteins. On the other hand, the other RFMs produce less proteins, as the pool that feeds these RFMs now has a smaller abundance of ribosomes. (3) We provide a detailed asymptotic analysis of the network steady state when the pool is starved. In particular, we show that in this case the initiation rates in every RFM become the bottleneck rates, and provide a closed-form expression for the total steady-state density and total production rate on any subset of RFMs in the network, relative to the steady-state pool density. (4) We provide an explicit bound for the total production rate of the network when the number of RFMs is very large. In particular, we show how the total density of ribosomes in the networks bounds the total production rate. This bound demonstrates that when the number of RFMs increases, the marginal utility of adding another RFM diminishes, and the total protein production saturates to a limiting value. (5) We demonstrate our analysis approach for the case of producing insulin proteins in a cell-free system, and show how it can provide useful guidelines for setting the parameters in such a system.
The remainder of this paper is organized as follows. The next section reviews the network model. Before going into the mathematical analysis, §3 illustrates several structural questions that can be studied using simulations of the model. Section 4 details the main mathematical results. Section 5 reports analyses of a biological system (gene expression in a cell-free system) with our models. For the sake of readability, all the proofs are placed in the appendix A. Throughout, we try to describe the biological implications of every analysis result.

The mathematical model
We use a model that includes m RFMs interconnected via a pool of free ribosomes. Each RFM models the dynamics of ribosome flow along one mRNA molecule and, in particular, each RFM may have a different length and different parameters (i.e. different codon decoding rates and initiation rates). The pool of free ribosomes represents ribosomes in the cell that are not attached to any mRNA. We begin by reviewing the various components of this network.

Ribosome flow model
The RFM includes n state variables x 1 , …, x n representing the normalized ribosome density in n sites along the mRNA, where each site corresponds to a group of consecutive codons. The density is normalized such that x i (t) ∈ [0, 1] for all t, where x i (t) = 0 represents that the site is empty, and x i (t) = 1 represents that the site is completely full. Thus, x i (t) may also be interpreted as the probability that site i is occupied at time t. The RFM also includes n + 1 positive parameters λ 0 , …, λ n , where λ i controls the transition rate from site i to site i + 1. In particular λ 0 controls the initiation rate, and λ n controls the exit rate.
The dynamics of the RFM is described by n balance equations, . . .

and _
x n ¼ l nÀ1 x nÀ1 ð1 À x n Þ À l n x n : To explain this, consider the equation for the change in density in the second site, namely, The term λ 1 x 1 (1 − x 2 ) represents the flow of ribosomes from site 1 to site 2. This is proportional to the transition rate λ 1 , to the density of ribosomes x 1 in site 1 and to the 'free space' (1 − x 2 ) in site 2. In particular, if site 2 fills up, i.e. x 2 is close to one, then the flow into site 2 decreases to zero. This is a 'soft' version of the simple exclusion principle, i.e. the notion that two particles cannot be in the same place at the same time. Similarly, the second term on the right-hand side of (2.2) is the flow from site 2 to site 3. Thus, equation (2.2) states that the change in density in site 2 is the flow from site 1 to site 2 minus the flow from site 2 to site 3. The exit rate from the last site is R(t) := λ n x n (t), and this is also the protein production rate at time t (see figure 1). Note that x i is dimensionless, and that λ i has units of 1/time. In all the biological simulations below, λ i is in units of 1/s. The state space of the RFM is the n-dimensional cube [0, 1] n . Since this invariant set is compact and convex, the RFM admits an equilibrium point e ∈ [0, 1] n . At an equilibrium, the flows into each site and out of each site are equal, so all the site densities are constant. Analysis of the equations describing the equilibrium point shows that there is a single equilibrium point e ∈ [0, 1] n (see [40]).
The RFM has been extensively used for studying the translation of a single, isolated mRNA. The model is highly amenable to analysis using tools from systems and control theory. It was shown in [40] that the RFM is a totally positive differential system [41] and this implies that any solution of the RFM converges to the unique equilibrium e. In particular, the protein production rate R(t) = λ n x n (t) converges to the steadystate production rate R := λ n e n . In other words, the positive transition rates λ 0 , …, λ n determine a unique steady-state density x 1 = e 1 , …, x n = e n along the mRNA, and for any initial density the dynamics converges to this profile.
Poker et al. [22] derived a useful spectral representation for the mapping from the rates λ 0 , …, λ n to the steady-state e. Given the RFM, consider the (n + 2) × (n + 2) tri-diagonal matrix :

ð2:3Þ
Since A is symmetric, all its eigenvalues are real. Since A is an irreducible matrix, with all entries non-negative, the Perron-Frobenius theorem [42] implies that A admits a simple maximal eigenvalue σ > 0, and the corresponding eigenvector z [ R nþ2 is unique (up to scaling) and satisfies ζ i > 0 for all i ∈ {1, …, n + 2}. Then, the entries of e satisfy [28] l n x n l n -1 x n -1 (1x n ) x 1 x 2 x n Figure 1. Ribosome flow model: arrows indicate the inflow and outflow into each site; x i is the ribosome density in site i (see (2.1)).
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220535 and the steady-state production rate satisfies In other words, the Perron eigenvalue and eigenvector of A provide all the information needed to determine the steadystate profile e, and the steady-state production rate R in the RFM. This spectral representation has several implications. For example, it implies that it is possible to numerically calculate efficiently the steady state even for very large RFMs using algorithms for computing the Perron eigenvalue and eigenvector of tri-diagonal matrices. Also, it follows from (2.5) that the function R = R(λ 0 , …, λ n ) is strictly concave [22], thus allowing to show that general steady-state protein production optimization problems are convex optimization problems [22]. It also reduces the sensitivity analysis of R with respect to any rate λ i to an eigenvalue sensitivity problem for the matrix A [21].

Ribosome flow model with input and output
As noted above, a cell typically includes a large number of mRNA molecules and ribosomes, that compete for shared resources, and studying translation on a single, isolated mRNA may thus provide limited insight on large-scale translation in the cell. To model translation in the cell requires a network of interconnected RFMs. The first step in building such a network is adding an input and output to the RFM. This yields the ribosome flow model with input and output (RFMIO) [34]. The RFMIO dynamics is ð1 À x n Þ À l n x n and y ¼ l n x n : ð2:6Þ The scalar input u : R þ ! R þ represents the density of ribosomes in the vicinity of the initiation site. Thus, if u(t) is large then the effective initiation rate at time t, given by u(t)λ 0 , increases. The scalar output y(t) = λ n x n (t) is the rate of ribosomes exiting the mRNA at time t. The additional input and output allow to connect RFMIOs in a network. Note that for u(t) ≡ 1, equation (2.6) reduces to the RFM.

The network
Raveh et al. [34] introduced a model composed of m RFMIOs interconnected via a pool of free ribosomes (see figure 2). Let n i , i = 1, …, m, denote the length of RFMIO #i. The state variables in RFMIO #i are denoted by x i 1 , . . . , x i n i . The density in the pool at time t is modelled by the scalar function z(t). The ribosomes that initiate translation in RFMIO #i are supplied from the pool through the pool output function G i (z(t)).
when the pool is empty the initiation rate in the RFMIO is zero), and G i (z) is continuous and strictly increasing in z (i.e. an increase in the pool density yields an increase in the initiation rates). Many possible functions satisfy these constraints, e.g. G i (z) = cz, with c > 0, and the uniformly bounded function G i (z) = αtanh (βz), with α, β > 0.
The pool feeds all the RFMIOs, and is fed by the ribosomes exiting all the RFMIOs, so the balance equation for the change in z(t) is where y i is the ribosome exit rate from RFMIO #i. Let denote the total density of ribosomes in the network at time t.
Since ribosomes cannot leave nor enter the network, In other words, s(t) is a first integral of the dynamics. It is important to note that there is no direct link between the RFMIOs in the network, and that the competition is not 'encoded' by changing the dynamical equations as was done in other context-aware models (see e.g. [43,44]). Competition arises only due to the interconnection via the shared pool of free ribosomes.
Raveh et al. [34] used results on cooperative systems with a first integral (see e.g. [37,38]) to prove that any solution of the network converges to a steady state. More precisely, for any p ≥ 0, let L p denote the p level set of the first integral, The network includes m RFMIOs interconnected via a pool of free ribosomes. For illustration only, we assume that RFMIO #1 has length n 1 = 2 and write its equations explicitly.
royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220535 that is, L p includes all the initial conditions of pool and RFMIO densities with total density s(0) = p. Then L p includes a single equilibrium point, and any trajectory emanating from an initial condition in L p converges to this equilibrium point. In other words, any two initial conditions of the network with the same total ribosome density will converge to the same equilibrium state. At this state, the total density is distributed along the mRNAs and the pool such that the flow into and out of each site is equal. Let e z ∈ [0, s(0)] denote the steady-state pool density, and let e i j denote the steady-state density in site j in RFMIO #i. Also, let i.e. e collects all the steady-state values in the network.
Here, we use the same network model as in [34] to study a different problem, namely, the effect of competition between multiple RFMIOs for scarce shared resources on the global dynamics of translation in the cell.

Simulation results
We begin with several synthetic simulation results that demonstrate the wealth of biological questions that can be addressed using the network model when allowing the number of RFMIOs m to vary, that is, when mRNAs are added or removed from the network. This also illustrates our general analysis approach that combines the spectral representation of the steady state in every RFMIO with the equation describing the first integral of the network.
We begin by considering a network that includes m identical RFMIOs, where each RFMIO has length n i = 2, i = 1, …, m. We also assume that every RFMIO is homogeneous, with λ 0 = λ 1 = λ 2 = 1. (Note, however, that all the theoretical results in §4 below hold for general lengths and rates.) We also assume that G i (z) = z for all i (i.e. the effective initiation rate is proportional to the number of free ribosomes in the pool).
To apply the spectral approach to each RFMIO in the network, let . Note that this is exactly the matrix (2.3) with n = 2, λ 0 G(e z ) = e z , and and the corresponding Perron eigenvector is It follows from (2.4) that the steady-state densities in every RFMIO are and since λ 2 = 1, the steady-state production rate of each RFMIO is RðcÞ ¼ e 2 ðcÞ: The equation for the total density of ribosomes s in the network is Combining this with (3.2) provides an explicit expression for s as a function of c. This can be inverted (at least numerically) to conclude for every total density s the corresponding c (and thus e z ), and then the spectral approach allows to obtain all the steady-state profiles in all the RFMIOs. The network allows to study how important steady-state quantities depend on the number of RFMIOs in the network. We first define several such quantities. The ratio between the density of ribosomes in the pool and the total density of ribosomes in the network is The total protein production at steady state, denoted TPR, is the production rate of all the RFMIOs in the network. Since here there are m identical RFMIOs, As m is increased, more ribosomes are attached to mRNAs and thus we can expect the steady-state pool density e z to go to zero. Then the initiation rate λ 0 G(e z ) = e z in each RFMIO becomes the bottleneck rate, and thus e i ≈ e z for i = 1, 2, in every RFMIO. Substituting this in (3.3) gives e z ≈ s/(1 + 2m), and the total production rate is then Thus, for a large m we expect the total steady-state production rate in the network to converge to s/2. Figure 3 depicts the exact network steady-state values for s = 50 as a function of the number of RFMIOs m. It may be seen that: (i) the steady-state densities e 1 and e 2 in every RFMIO decrease monotonically with m. The same holds for q = e z /s, and (ii) the total production rate me 2 increases monotonically with m, and converges to a saturation value of s/2 = 25.
Summarizing, the addition of mRNAs increases the total protein production, as ribosomes now also translate the new mRNAs. However, it decreases the translation rate of the other mRNAs by depleting the pool of free ribosomes. As the number of mRNAs increases, the marginal utility in adding another mRNA decreases to zero. The total production rate in the network is bounded, and an important factor in this bound is the total density of ribosomes in the network.
These simulation results suggest that when optimizing the production rate of a synthetic system (e.g. in highthroughput heterologous gene expression) it may be useful to use mRNA levels below the regime where the marginal utility in adding another mRNA become negligible.
The next section provides a rigorous analysis of these topics.

Main results
We begin by considering the effect of adding an RFMIO to the network. We introduce some notation. Recall that RFMIO #i is characterized by the tuple where n i is the length of RFMIO #i, G i : R þ ! R þ is the ith pool output function, and l i j are the rates along RFMIO #i. Consider a network with m + 1 RFMIOs obtained by adding an RFMIO to a network of m RFMIOs. Let e z (m) [e z (m + 1)] denote the pool density in the network with m [m + 1] RFMIOs. What is the relation between the steadystate pool densities before and after adding RFMIO #(m + 1)? The next result shows that e z (m + 1) is always smaller than e z (m) (cf. figure 3a). The proof is placed in appendix A.1. In other words, adding an RFMIO to the network, while keeping the total ribosome density constant, always decreases the steady-state pool density. This makes sense, as the new mRNA 'consumes' ribosomes from the pool. exists. The next result shows that this limit is zero. Since we take m → ∞, we need to impose some technical conditions on the RFMIOs.
Assumption 4.2. From here on we always assume that the following properties hold.
(1) There exists l Ã . 0 such that i.e. all the initiation rates are bounded from below by l Ã . (2) There exists λ* > 0 such that i.e. all the exit rates are bounded from above by λ*.
(3) There exist p > 0 and g Ã . 0 such that for any z ∈ [0, p], we have i.e. all the pool output functions G i (z) are bounded from below by the linear function g Ã z on the interval [0, p].
These three conditions are clearly reasonable. The next result shows that under these conditions, the limit in (4.2) is zero (cf. figure 3a). The proof is placed in appendix A.2. From a biological point of view, the case m → ∞ may seem unreasonable. However, Proposition 4.3 implies that for a large m, e z (m) will be small. This scenario is the focus of this paper. Indeed, when the pool includes a large number of free ribosomes there is little competition, and thus the indirect coupling between the mRNAs is weak. The interesting case is thus when the pool is close to being depleted, that is, when e z is small. In the biological context, this may represent a cell under stress conditions, or under a high-yield viral infection, or when there are many ribosomal 'traffic jams' along the mRNAs, so the pool is depleted.
The next result uses the spectral representation to analyse the steady-state densities and production rate in every RFMIO when the pool is starved. royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220535 Proposition 4.4. Consider a network of RFMIOs T 1 , T 2 , . . ., interconnected via a pool. Suppose that either the total ribosome density s goes to zero, or that the number of RFMIOs tends to infinity, so that e z → 0. Then for any i, RFMIO #i satisfies In particular, if the pool output function G i is differentiable at zero, then The proof is placed in appendix A.3.
In other words, when the pool becomes depleted (e.g. because m is large or the total ribosome density s is small) every density along the ith mRNA behaves asymptotically like the pool output function G i (e z ). This makes sense, as the effective initiation rate l i 0 G i ðe z Þ becomes the bottleneck rate in the mRNA. Note that (4.4) implies that e i j =l i 0 G i ðe z Þ is inversely proportional to l i j . This is reasonable, as l i j controls the flow out of site j. These values are taken from Zarai et al. [45] and correspond to the S. cerevisiae gene YBL025W that encodes the protein RRN10, which is related to regulation of RNA polymerase I. This gene has 145 codons (excluding the stop codon), and was divided into six consecutive groups of codons: the first group includes the first 24 codons (that are also related to later stages of initiation). The other groups include 25 nonoverlapping codons each, except for the last one that includes 21 codons. This partitioning was found to optimize the correlation between the RFM prediction and biological data measurements (see [45] for more details). We increased m from 1 to 1000 while keeping the total ribosome density fixed at s = 20, thus depleting the steady-state pool density. The input functions are G i (x) = x for all i. Since all the RFMIOs are identical, it is sufficient to consider the steadystate density in a single RFMIO. Figure 4 depicts e j /e z , j = 1, …, 5, as a function of (1/e z ). It may be seen that as e z decreases, every ratio e j /e z converges to the asymptotic value given in proposition 4.4.
The following result is an immediate corollary of proposition 4.4. Recall that the constant total ribosome density in the network is and the total steady-state production rate is Also, let q := e z /s denote the ratio between the free ribosomes and the total number of ribosomes in the network. (2) The ratio between the free ribosomes to the total number of ribosomes in the network satisfies where Hðl i 1 , . . . , l i n i Þ : ¼ n i ð P n i j¼1 ð1=l i j ÞÞ À1 is the harmonic mean of the rates l i 1 , . . . , l i np , that is, all the rates except for the initiation rate.
The proof is placed in appendix A.4. These results provide closed-form asymptotic expressions for important biological quantities when the pool is starved. Note that as e z → 0, TPR/e z depends on all the initiation rates l i 0 , but not on any of the other rates (since the number of ribosomes along any mRNA is low, there are no traffic jams). However, the ratio between the density of ribosomes in the pool and the total number of ribosomes in the network does depend on the harmonic mean of all the rates in the network.
We note that a similar closed-form expression can be obtained for the total steady-state density and total production rate on any subset of RFMIOs in the network.   royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220535 Corollary 4.6 implies that when the pool is starved we can replace an entire network of m identical RFMIOs by a much simpler network, while keeping the asymptotic steady-state properties unchanged. Proposition 4.8. Let s > 0 denote the total number of ribosomes. Consider the following two networks of RFMIOs: (1) A network of m identical RFMIOs, each of length n, rates λ 0 , …, λ n , and output functions G(z) = gz, with g > 0. Let TPR [q] denote the total production rate [ratio between free ribosomes and s] in this network. (2) A network consisting of a single RFMIO of length n ¼ 1, rates l 0 , l 1 , and GðzÞ ¼ gz, with g . 0. Let TPR [q] denote the total production rate [ratio between free ribosomes and s] in this network.
If the parameters of the second network are chosen such that The proof is placed in appendix A.5.
In other words, we can use a network with a pool and a single RFMIO, with a single site, to simulate and analyse the first network. Note that conditions (4.10) and (4.11) are quite intuitive. Roughly speaking, the first implies that the initiation rates in the two networks are equal (taking into account that in the first network there are m RFMIOs and in the second a single RFMIO), whereas the second condition requires that some mean of the other rates along the RFMIO is also equal. Recall that this has m identical RFMIOs of length n = 5 and the rates given in (4.6). In this case, g = 1, λ 0 = 0.1678 and X 5 j¼1 ðl j Þ À1 ¼ 18:6512: Proposition 4.8 implies that we can replace this network of m RFMIOs by a network consisting of a single RFMIO of length one, with GðxÞ ¼ x, l 0 ¼ 0:1678 m and l 1 ¼ 1=18:6512, and the asymptotic behaviour of the two networks when the pool is starved will be identical.

A biological example: production of insulin in a cell-free system
In this section, we apply our model to analyse the production of insulin protein in a cell-free system. A cell-free system is an in vitro-based approach to study and/or generate biological reactions that take place within cells via the isolation of relevant cellular components (e.g. ribosomes, RNA polymerases, tRNA molecules). This approach reduces the complex interactions typically found when working with whole cells. Cell-free systems are often used in biotechnology for heterologous gene expression [46]. Data on the coding region of insulin fitted to S. cerevisiae was taken from Kjeldsen [47]. The codon decoding rates, that are based on the analysis of typical decoding rates in vivo from ribo-seq data [48], were taken from Dana & Tuller [49]. Groups of 10 consecutive codons were coarse grained into one RFMIO site, as was done in previous studies (see e.g. [31]). This yields 12 sites. The transition rate λ i of site i is the inverse of the sum of the decoding times along the 10 related codons. These rates were then normalized such that their average value is 10 codons per second (the typical decoding rate in eukaryotes [50]). 1 Since we study the translation of insulin in a cell-free system, we assume that all the mRNAs in the system encode the protein insulin. Thus, the network includes m identical RFMIOs interconnected via a pool.
We set s = 25 × 10 4 , m = 39.5 × 10 3 (see the data in [31,51,52]) and assume that all the pool output functions are identical: G i (x) = cx, for all i = 1, …, m. To calibrate the constant c, we assume that c maximizes the effective initiation rate while keeping it lower than all the other transition rates. Mathematically, this yields the equation and this gives c = 7.4758 × 10 −5 .
To study the scenario where the pool is starved, we consider two cases: in the first we fix m and decrease s, and in the second we fix s and increase m. From a biological perspective, both cases correspond to the fact that ribosomes may be 'more expensive' then mRNAs, and the goal is to optimize production while using a minimal number of ribosomes.

Varying the total density of ribosomes
Consider the case where m = 39.5 × 10 3 is fixed, and s varies. Figure 5 depicts TPR/e z (that is, the ratio between the steadystate total production rate and the steady-state pool density) as a function of m/s (that is, the number of mRNA molecules divided by the total number of ribosomes in the network). As

Varying the number of mRNAs
Consider the case where s = 25 × 10 4 is fixed, and m is varied. Figure 6 depicts the number of free ribosomes e z as a function of m/s. As m increases, e z sharply decreases, as more ribosomes attach to the additional mRNAs, and thus the pool is starved. In particular, e z goes to zero as m → ∞ (see proposition 4.3). In practice, when m/s = 80 we already get a very low value of e z , and then the asymptotic results described in our analysis can already be used. Note that these results can provide important guidelines for setting the cell-free system parameters. For example, figure 5 shows that to achieve a steady-state production rate that is 99% of the maximal possible production rate we should set m/s = 10, that is, 10 mRNA molecules for each ribosome in the system. Our model allows to estimate the total production rate for any m/s value.

Discussion
The competition for shared resources plays an important role in gene expression. It is known that competition for free RNAPs and free ribosomal subunits is a major bottleneck for gene expression in bacteria [53]. Competition for shared resources also hampers our ability to reliably predict the behaviour of synthetic biology constructs (see [13,43,54,55] and references therein).
We considered a network of mRNAs fed by a pool of free ribosomes in the scenario when the pool is starved. This scenario is expected to be relevant for example in biological networks under stress conditions or viral attack, and in synthetic networks designed to optimize the total production rate.
We used a mathematical model of a network of RFMIOs connected via a pool of free ribosomes. Using the spectral representation of the RFMIO steady state we derived closed-form expressions for several relevant biological quantities in the regime when the pool is starved. These include the total protein production rate in the network, and the ratio between the number of ribosomes in the pool and the total number of ribosomes in the network.
We demonstrated the analytical results using both synthetic examples and an example based on biological data. The results reported here can be used both in systems biology studies of natural systems and for the design of synthetic constructs. We believe that an interesting direction for further research is to use our results in the biological context of a cell attacked by viruses. This will allow to study both qualitatively and quantitatively important questions. For example, can the virus effectively shut down the host protein production (and in particular immune-related proteins) by depleting the pool of free ribosomes or are other mechanisms needed?
Currently, there are no systematic measurements of the dynamic ratio between the number of free ribosomes and the number of mRNAs in cells. Nevertheless, stress conditions are known to induce ribosomal traffic jams and thus we may expect a significant reduction in the number of free ribosomes. Even in the presence of physiological feedback on the number of ribosomes and mRNAs in such conditions, production of new ribosomes may be slow as it consumes considerable cellular resources. Novel experimental procedures are needed to study these issues in the cell under various conditions.

Appendix A. Proofs
For the sake of readability, we begin with a short and general description of our analysis approach. Consider without loss of generality the steady-state densities e 1 ¼ ½e 1 1 . . . e 1 n1 T along RFMIO #1 in the network. The spectral representation implies that we can retrieve e 1 from the Perron root and Perron eigenvector of the matrix royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220535 Indeed, at steady state the pool density is z(t) = e z , so the initiation rate at RFMIO #1 is l 1 0 G 1 ðzðtÞÞ ¼ l 1 0 G 1 ðe z Þ. When e z is close to zero, so is G 1 (e z ) and this implies that entries (1, 2) and (2, 1) in A 1 are very large. We use an asymptotic analysis of the spectral properties of A 1 to derive an approximate expression for e 1 and, similarly, for any e i , i = 1, …, m. Now by (2.8) and (2.9), and thus we obtain the entire steady state of the network. We begin with several auxiliary results that describe the asymptotic spectral properties of a specific tri-diagonal matrix. We use R n þ : ¼ fx [ R n : x i ! 0, i ¼ 1, 2, . . . , ng to denote the non-negative orthant in R n , and R n þþ : ¼ fx [ R n : x i . 0, i ¼ 1, 2, . . . , ng to denote the positive orthant in R n . For a Hermitian matrix S ∈ C N×N we denote its eigenvalues by Recall also that if j Á j : C N ! R þ is a vector norm, and k Á k : C NÂN ! R þ is the induced matrix norm, then |σ i-(S)| ≤ kSk for any i.
We can now prove theorem A.1. First note that for any c ≥ where the last inequality follows from the fact that This completes the proof of (A 2). Taking c → ∞ in (A 2) proves (A 3).

c). This proves (A 5). ▪
In practice e z is positive (even if small), so it is useful to derive explicit error bounds in the asymptotic expressions. Note that since we are interested in ratios between entries of the Perron vector, we do not necessarily require it to be normalized. For two functions f, g : R þ ! R þ , we write f = O(g) if there exist c, y > 0 such that jf ðxÞj cjgðxÞj for any x ! y: ðA 15Þ The next result provides more explicit information on the Perron eigenvector of T[c, α].  Note that this implies that R does not depend on c, and that Let ν(c) be the Perron eigenvector of S(c) defined in (A 10), and let Pð1Þ : ¼ Proj sp(nð1Þ) and PðcÞ : ¼ Proj sp (nðcÞ) ðA 19Þ be the projection operators on sp(ν(∞)) and sp(ν(c)), respectively. Let jðcÞ : ¼ PðcÞnð1Þ. Since we project on sp(ν(c)), this implies that ξ(c) is a Perron eigenvector of S(c) corresponding to σ(c)/c. Let η(c) := D −1 ξ(c). Then (A 11) implies that η(c) is a Perron eigenvector of T(c) corresponding to σ(c). We will show that η(c) satisfies (A 16). We begin by analysing ξ(c). Let G be the circle in the complex plane parametrized by gðtÞ ¼ 1 þ 1 2 expðitÞ, t ∈ [0, 2π) (see figure 7).
Using Weyl's inequalities (A 6), the second largest eigenvalue of T(c) satisfies σ 2 (T(c)) ≤ 2kRk. The proposition then follows from (A 12 We require the following auxiliary result that provides a Neumann series representation for the difference between these two projections. The proof is provided for completeness. Let i.e. the maximal norm of the resolvent of S(∞) over G. This maximum exists since (λI N − S(∞)) −1 is continuous in a neighbourhood of G.
In particular, the series on the right-hand side of (A 21) converges absolutely and uniformly on G, and royalsocietypublishing.org/journal/rsif J. R. Soc. Interface 19: 20220535 ðc À2 RðlI N À Sð1ÞÞ À1 Þ k , and the series converges if kc −2 R(λI N − S(∞)) −1 k < 1. We have kc À2 RðlI N À Sð1ÞÞ À1 k qðcÞ : ¼ c À2 M G kRk, ðA 23Þ and q(c) < 1 for any c . maxð ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi M G kRk p , 4kRkÞ. Moreover, for any such c, we have  Remark A.5. Note that approximations of the coordinates of η(c) up to arbitrarily small error can be obtained by computing higher order terms in the Neumann series (A 22). As an example, consider for simplicity the case α i = 1 for i = 1, …, N − 2, which implies kRk = 1. Assume that c . maxð ffiffiffiffiffiffiffi M G p , 4Þ and let δ > 0 be the desired error bound. Let L > 0 be a sufficiently large integer such that   can be computed explicitly by the Cauchy residue theorem [57]. Using the fact that η(c) = D −1 ξ(c) yields so the entries of D −1 μ(c) approximate those of η(c) with an error smaller than δ.
The next result provides an explicit expression for (λI N − S(∞)) −1 with l [ G. This can be used to derive an upper bound on the constant M G , given in (A 23), thereby leading to an explicit estimate of the errors in (A 29). Moreover, (λI N − S(∞)) −1 can be substituted into (A 31) to obtain highorder approximations of the Perron eigenvector.
Proposition A.6. Let G be as in figure 7. Then, for any l [ G and B defined in (A 33). Since λ ≠ 0, it can be easily verified that C −1 = E and A −1 = D. Hence, and this completes the proof. ▪ We can now prove the results in §4.
statement of the proposition and (A 43), and the fourth from the definition of 1. However, since the system is initialized at the steady state, P n p j¼1 _ x p j ¼ 0 for all t ≥ 0. This contradiction completes the proof.

A.4. Proof of corollary 4.6
Assume a sequence of networks with m p RFMIOs, p = 1, 2, …, with lim p→∞ m p = ∞. For any p, let e z ( p) be the corresponding steady-state pool density. By proposition 4.3, lim p→∞ e z ( p) = 0. For any p, let f p ðxÞ : ¼ P mp i¼1 ðl i n i e i n i ðpÞ=e z ðpÞÞ1 ½i,iþ1Þ ðxÞ. Note that f p (x) ≥ 0 for all x, and that Ð 1 À1 f p ðxÞ dx ¼ P mp i¼1 ðl i ni e i ni ðpÞ=e z ðpÞÞ is the ratio between the total production rate and the pool density for the case of m p RFMIOs in the network. By proposition 4.4, the functions ff p g 1 p¼1 converge pointwise to f ðxÞ : The last equality follows from assumption 4.2, which implies that l i 0 G 0 i ð0Þ ! l Ã g Ã . 0. In the case where the number of RFMIOs is fixed at m as e z → 0, by replacing Fatou's lemma with the dominated convergence theorem [58, Theorem 2.8.1], the limit-inferior can be replaced by a limit, and the inequality can be replaced by equality. Note that in this case P m i¼1 l i 0 G 0 i ð0Þ , 1. This proves (4.8). To prove (4.9), consider first the case where the number of RFMIOs m is fixed and the total density in the network satisfies s → 0. Recall that q = e z /s, and (4.7) gives q ¼ ð1 þ P m i¼1 P n i j¼1 ðe i j =e z Þ À1 . Using , and this completes the proof. We note in passing that in the case where the total ribosome density in the network s is fixed and the number of RFMIOs satisfies m → ∞, assumption 4.2 gives A.5. Proof of proposition 4.8 Using (4.8) gives lim ez!0 ðTPR=e z Þ ¼ l 0 gm and lim ez!0 ðTPR=e z Þ ¼ l 0 g, and combining this with (4.10) gives (4.12). Similarly, equation (4.9) gives lim ez!0 q ¼ gl 0 mn Hðl 1 , . . . , l n Þ and lim ez!0 q ¼ gl 0 n Hðl 1 , . . . , l n Þ , and using (4.10) and (4.11) implies that these expressions are equal. This completes the proof of proposition 4.8.