Mediating Ribosomal Competition by Splitting Pools
MMediating Ribosomal Competition by Splitting Pools
Jared Miller, M. Ali Al-Radhawi, and Eduardo D. Sontag
Abstract
Synthetic biology constructs often rely upon the introduction of “circuit” genes into host cells, in orderto express novel proteins and thus endow the host with a desired behavior. The expression of these newgenes “consumes” existing resources in the cell, such as ATP, RNA polymerase, amino acids, and ribosomes.Ribosomal competition among strands of mRNA may be described by a system of nonlinear ODEs called theRibosomal Flow Model (RFM). The competition for resources between host and circuit genes can be amelioratedby splitting the ribosome pool by use of orthogonal ribosomes, where the circuit genes are exclusively translatedby mutated ribosomes. In this work, the RFM system is extended to include orthogonal ribosome competition.This Orthogonal Ribosomal Flow Model (ORFM) is proven to be stable through the use of Robust LyapunovFunctions. The optimization problem of maximizing the weighted protein translation rate by adjusting allocationof ribosomal species is formulated and implemented.
I. I
NTRODUCTION
The process of protein expression is mediated by ribosomes. After a gene in DNA has been transcribedinto messenger RNA (mRNA), a ribosome binds to the mRNA and begins protein translation. ThemRNA is divided into a set of 3-nucleotide segments called codons, and each codon corresponds to anamino acid or a stop instruction. The ribosome attracts a tRNA carrying an amino acid which matchesthe currently read codon, and appends it to a growing polypeptide chain. Once the ribosome hits a stopcodon, it falls off the mRNA and releases the amino acid chain as a polypeptide, which is subsequentlypost-translationally processed into a final protein product. The mRNAs remain in the cell until theyare degraded or destroyed (such as by siRNA or RISC). Intact mRNAs continually attract ribosomesand produce protein, and all mRNAs in the cell compete for a finite pool of resources which includesribosomes.One way to decouple host and circuit genes is to split the pool of resources via orthogonal ribosomes[1], [2]. Orthogonal ribosomes are mutated ribosomes which only translate specifically modified circuitgenes, and can be constructed by introducing a synthetic 16s rRNA into the cell. The circuit’s mRNAsbinding sites are designed so that only mutated ribosomes will translate circuit mRNAs; host ribosomesare not attracted to the circuit mRNAs [3]. Orthogonal ribosomes may be able to decrease competitionand increase protein throughput by splitting the pool.We present the Orthogonal Ribosomal Flow Model (ORFM) extending the existing RFM. The ORFMsystem is a conservative system where the total number of ribosomes is conserved. Global asymptoticstability of the ORFM is certified by means of robust Lyapunov functions. A simple bisection algorithmis detailed to compute the unique ORFM steady state. The protein throughput of the ORFM systemcan be changed by adjusting the production rate of different species of ribosomes, and maximizing thisthroughput is a non-concave optimization problem.The structure of the paper is as follows: Section II reviews existing models for protein translationand ribosome competition. Section III introduces the ORFM. Section IV adds a feedback controller toregulate production of ribosomes. Section V presents a non-concave optimization problem to maximizeprotein throughput. Section VI concludes the paper. The stability of the ORFM is proved in the Appendix.
The authors are with the Department of Electrical & Computer Engineering, Northeastern University, Boston, MA. E.D. Sontag, is alsowith the Department of Bioengineering, Northeastern University, and the Laboratory of Systems Pharmacology, Harvard Medical School,Boston, MA. Emails: { miller.jare,malirdwi,e.sontag } @northeastern.edu . a r X i v : . [ q - b i o . M N ] S e p ig. 1: A pictorial representation of the inflow, outflow, and densities of a 5-codon RFMIO at times t = 1 . and ∞ . II. R IBOSOMAL F LOW M ODELS (RFM S ) A. Single-Strand Translation Models
The first models analyzing steady-state ribosomal distributions on mRNA were based on lattice models.In a Totally Asymmetric Simple Exclusion Process (TASEP), mRNA is represented as a finite-lengthone-dimensional lattice of codons [4], [5]. Each ribosome waits a random amount of time (exponentiallydistributed) before attempting to move rightwards to the next codon.The RFM is a deterministic mean-field approximation to the TASEP model [6]. Each codon on themRNA has a normalized ribosomal density x j ( t ) ∈ [0 , , which one may think of as the probabilitythat a ribosome is present on codon j at time t . Transition rates between codons are denoted here as λ j .The initiation rate λ is the rate at which the mRNA attracts the ribosome to begin translation, whilethe λ j
Rates λ , Constant Input u output: Codon Steady States e j µ j = 1 / (cid:112) λ j for j = 0 . . . nJ = tridiag ( n +2 , µ ) σ = σ max ( J ) , ζ = v max ( J ) e j = µ j σ ζ j +2 ζ j +1 for j = 1 . . . n The RFM system is a monotone control system, and the set of admissible controls u ∈ U is the set ofbounded and measurable functions u ∈ R + . The RFM is also state and output-controllable, and desiredtranslation rates and patterns can be achieved by proper choice of u and λ [8]. B. Ribosomal Competition
In real genetic systems, protein translation on a strand of mRNA does not occur in isolation. At anyparticular moment, all mRNA’s in the cell compete for a finite (and possibly time-varying) numberof ribosomes. Particularly slow transition rates λ on codons can lead to strands of mRNA hoardingribosomes, reducing the availability of ribosomes in the pool for all other mRNA’s and leading to aglobally depressed translation rate.Competition was introduced into the RFM framework in the case of one strand of mRNA by positivefeedback [9] and on circular mRNA [10], in each case resulting in systems with a globally asymptoticallystable steady state. Raveh et. al. introduced a ribosomal flow model network with a pool (RFMNP)to abstractly describe the impact of ribosomal competition [11]. Each of the s strands of mRNA ismodeled by a RFMIO ( x ij for codon j of mRNA i ), and all RFMIO are connected to a common pool z . The total number of ribosomes in the system N r = z ( t ) + (cid:80) i,j x ij ( t ) is a conserved quantity. Theinput u of each mRNA is a monotonically increasing concave saturation function G i ( z ) (commonly z or tanh ( z ) ), which describes the likelihood that a ribosome from the pool will attach itself to strand i . The output y i = λ in x in is the translation rate, and ribosomes leaving x i return to the pool. The pooldynamics are therefore: ˙ z = (cid:80) i λ in x in − (cid:80) i λ i G ( z )(1 − x i ) (2)The RFMNP is a closed loop system. The number of ribosomes N r defines a stoichiometric class,and RFMNP has a globally asymptotic equilbrium point with all e ij ∈ (0 , and z ∈ (0 , N r ) foreach N r . Stability can be proven by contraction, where the weighted L norm between trajectories isnon-expanding over time [11].Competition effects can be observed by perturbing parameters λ ij and analyzing resultant steady states.If λ ij changes for a specific strand of mRNA x i , then the protein translation rate y i will change in thatsame direction (increasing λ ij will increase y j ). The translation rate of all other mRNA in the networkwill change uniformly: an increase in λ ij will raise y i and will either raise or lower all y k (cid:54) = i . Furtherdiscussion on competition effects can be found in Section 4.2 of [11], and some of these results aredemonstrated in Figure 2. Figure 2 shows four examples of an RFMNP where each of the three strandsof mRNA have five codons each. All strands start with homogeneous transition rates of λ ij = 5 , and thereare a total of N r = 5 ribosomes in the system. Figure 2a visualizes the steady state of this homogenoussystem, where each mRNA has a protein translation rate of y ss = 1 . . Ribosomes cannot easily passfrom codons 4 to 5 because λ = 0 . . This bottleneck drops the translation efficiency of strand 1 to y ss = 0 . and the others to y ss = 0 . . Figure 2c has λ = 40 , so strand 1 attracts ribosomes fromthe pool at a much faster rate than the other strands. The first codon cannot drain ribosomes quickly a) Homogeneous transitions. (b) Strand 1 is slow(c) Strand 1 is greedy (d) Strand 1 is fast Fig. 2: RFM competition at steady statesince λ ... = 5 , so ribosomes are stuck in strand 1. Translation rates rise for strand 1 to y ss = 1 . , andfall on other strands to y ss = 0 . . Figure 2d modifies strand 1 to λ = 40 . Ribosomes quickly exitstrand 1 and enter the pool again for use in further translation. This improves the translation efficiencyof all strands, raising y ss = 1 . and y ss = 1 . .III. O RTHOGONAL R IBOSOMAL F LOW M ODEL
Orthogonal ribosomes can be added to the RFM scheme (ORFM) by splitting the ribosome pool z . Codefor simulation and visualization is publicly available online at https://gitlab.com/jarmill/Ribosomes . A. ORFM Formulation
The proposed model of ribosome translation has M species of ribosomes. Each ribosome species hasa pool z p for p = 1 . . . M of available ribosomes. Strands of mRNA that use ribosome type p form anRFNMP with pool z p .Ribosomes of type p are formed by the combination of a os-15 rRNA of type p and the remainingribosome components. The protein backbone and large ribosomal subunit are treated as an ‘Empty’ribosome E , which has no translation capacity on its own. Empty ribosomes bind with rRNA at rate k + p , and the ribosomal complex dissociates at rate k − p , assuming an effectively infinite amount of rRNAavailable. These kinetics are inspired by the work of Darlington et. al. in their mass-action model ofprotein translation and cell metabolism [12]. The M pools of translating ribosomes z p are each connectedto the pool of empty ribosomes z E .Each of the s p strands of mRNA that obtain ribosomes from pool z p will have a corresponding length n pi and are repeated with multiplicity m pi . The codon at location j ∈ . . . n pi − of RNA strand ∈ . . . s p coming from pool p is x pij . The constants of this system are the N r ribosomes, binding rates k + p and k − p , and transition rates λ pij . The total number of ribosomes N r is a conserved quantity: N r = z E + (cid:80) Mp =1 z p + (cid:80) Mp =1 (cid:80) s p i =1 (cid:16) m pi (cid:80) n pi j =1 x pij (cid:17) (3)The translation dynamics are: ˙ z E = (cid:80) p k − p z p − (cid:80) p k + p z E (4a) ˙ z p = k + p z E − k − p z p + (cid:80) i m pi λ pin pi x pin pi − (cid:80) i m pi λ pi (1 − x pi ) G pi ( z p ) (4b) ˙ x pi = λ pi (1 − x pi ) G pi ( z p ) − λ pi (1 − x pi ) x pi (4c) ˙ x pij = λ pij − (1 − x pij ) x pij − − λ pij (1 − x pij +1 ) x pij (4d) ˙ x pin pi = λ pin pi − (1 − x pin pi − ) x pin pi − λ pin pi x pin pi (4e)Define ¯ N r = z E + (cid:80) p z p . For a fixed ¯ N r , the number of ribosomes in each pool can be obtained bysolving a linear system where K p = k + p /k − p : z p = K p z E ¯ N r = z E + (cid:80) p z p (5)Algorithm 2 computes the steady-state of the general model by bisection on ¯ N . This algorithm can Algorithm 2:
ORFM Steady State input : λ pij , N r , G pi , K p , (cid:15) output: Steady States z E , z p , e pij ¯ N min = 0 , ¯ N max = N r repeat ¯ N curr = ( ¯ N max + ¯ N min ) / z E = ¯ N r (cid:80) p K p z p = K p ¯ N r (cid:80) p K p e pi = RFMIO SS ( λ pi , G pi ( z p )) (Alg. 1) N curr = ¯ N curr + (cid:80) p,i,j e pij if N curr ≤ N r then z min ← z curr else z max ← z curr enduntil | N curr − N r | ≤ (cid:15) ;compute the steady state of an RFMNP by treating the network as a -species model with K → inf where ¯ N = z .Figure 3 shows examples of mRNA competing for N r = 10 ribosomes, where each mRNA may takeribosomes from either the host pool ( p = 1 , blue) or the circuit pool ( p = 2 , red). The ‘empty’ ribosomesare displayed in purple in the third tank. In this network, k + p = k − p = 1 , so at steady state, the poolquantities z = z = z E . In figure 3a no mRNA takes ribosomes from the circuit pool z , so theribosomes in z induce deadweight loss for the system.Figure 4 illustrates a general model orthogonal ribosomal system with pool across three ribosomalsubspecies and N r = 20. Now k + p = k − p = 0 . , p = 1 , , , so all pools will have an equal number ofribosomes at equilibrium ( z =0 . . a) No circuit demand (b) High circuit demand Fig. 3: Orthogonal RFM visualizationsFig. 4: ORFM with 3 ribosomal subspecies
B. Stability of the ORFM
The RFM has recently been studied by constructing explicit Robust Lyapunov Functions (RLFs) [13]via writing it as a Chemical Reaction Network (CRN) and utilizing relevant methods [14], [15], [16].Such techniques provide explicit formulae of Lyapunov functions for general kinetics (not limited toMass-Action kinetics), and have an easy-to-use software package.In this subsection, we derive the stability of the ORFM model (4a)-(4e) via an RLF. Let x =: [ x , , .., x M,s M n M,sM ] ∈ [0 , N c be the vector of all codon occupancies, and N c be the total number of codons. Let z :=[ z , .., z M , z E ] T ∈ [0 , N r ] M +1 be the vector of all pool occupancies. We therefore have: Theorem 1:
Consider the system (4a)-(4e). Then,1) The function: V ( x, z ) = (cid:88) p,i m pi (cid:88) j | ˙ x pij | + (cid:88) p | ˙ z p | + | ˙ z E | , is a (non-strict) Lyapunov function for any choice of { λ p,ij } > and monotone functions G pi , and(2) For any fixed total number of ribosomes in the system N r > , there exists a unique positive globallyasymptotically stable steady-state ( x e , z e ) for (4a)-(4e).A proof can be achieved by transforming (4a)-(4e) into a CRN and extracting a stoichiometry matrix.For a fixed number of pools, strands, and codons, Theorem 1 can be verified via the software package LEARN [13]. A proof for the general result is given in the Appendix.
C. Cross-talk Model
The analysis of the previous sections assume that only host ribosomes translate host genes and onlycircuit ribosomes translate circuit genes. While it is energetically unfavorable for host ribosomes totranslate circuit genes, this event may still happen at a low proability. The phenomenon of a ribosomeof species p translating an mRNA that prefers ribosomes of species p (cid:48) (cid:54) = p is a ’cross-talk’. Each codonof mRNA still can only accomodate one ribosome at a time, so ribosomes of different species areadditionally competing for space. x pj is the probability that a ribosome of species p sits at codon j ,and x j = (cid:80) p x pj is the probability that some ribosome is present at codon j . A cross-talk RFM can betreated as a MISO polynomial (non-tridiagonal) system, with inputs u p and output y = (cid:80) Mp (cid:48) =1 λ pn x pn : ˙ x p = λ p u p (1 − M (cid:88) p (cid:48) =1 x p (cid:48) i ) − λ p (1 − M (cid:88) p (cid:48) =1 x p ) x p ˙ x pj = λ pj − (1 − M (cid:88) p (cid:48) =1 x p (cid:48) ij − ) x pj − λ pj (1 − M (cid:88) p (cid:48) =1 x pj ) x pj +1 ˙ x pn p = λ pn p − (1 − M (cid:88) p (cid:48) =1 x pn − ) x pn p − λ pn x pn (6)Fig. 5: mRNA with 3 ribosomal speciesThere is no requirement that the different transition rates λ pij be the same between species. Figure 5shows a 5-codon mRNA with λ pj uniformly drawn from the range [0 . , . .A cross-talk RFM where of length n where λ pj = λ j , j (cid:54) = 0 is equivalent to a standard RFM in termsof sums s j = (cid:80) Mp =1 x pj , with effective initial rate λ s = (cid:80) Mp =1 λ pi . Ribosome flux in codons can beequivalently expressed as:ig. 6: A network of cross-talk mRNA ˙ x pj = λ j − (1 − s j ) x pj − − λ j (1 − s j +1 ) x pj (7)Equilibrium values for s can be found through spectral methods, and the steady state per-species codonvalues are: x pj In an open loop crosstalk model, the various z p are constant or are external inputs to the system. Thecrosstalk model can also have network competition, in which the empty-tank pool dynamics are added.As before, x pij is the probability that a ribosome of species p is located on codon j of mRNA i at anyparticular time. The input u pi applied to each strand is G pi ( z p ) , and pool dynamics ˙ z p are adjusted todeal with the cross-talk effects (1-sum occupancy).Figure 6 illustrates a cross-talk orthogonal RFM with N r = 10 ribosomes.IV. S ELF -I NHIBITING F EEDBACK C ONTROLLER This section introduces a feedback controller inspired by [12] to dynamically generate orthogonalribosomes as needed in an ORFM. Figure 3a shows an example where there are no mRNAs thattake ribosomes from the circuit pool z . A high k +2 therefore leads to the creation and maintenanceof circuit ribosomes that are not used in translation. A feedback controller can be introduced to createorthogonal ribosomes only as needed, and decrease wastage in the system. This can be accomplished bycreating one distinguished mRNA x pF that translates a protein F p for each non-Host pool. The dynamicsof the protein F p with a degradation rate δ p is: ˙ F p = λ pFn x pFn − δ p F p (9)A Hill-like inhibition term can be used to suppress the creation of new ribosomes (based on Equation 18of Supplement 1 of [12]). For a nominal ribosome creation rate k +0 p , exponent γ p , and constant F p , thesuppressed k + p is: k + p = k +0 p / (1 + F p /F p ) γ p . If δ = 0 , F p never degrades once translated. Asymptotically, a) No circuit demand (b) High circuit demand Fig. 7: ORFM with inhibitor protein F p . k + p → and z p → . A choice of δ > will lead to z p > , but proper choice in inhibitor parameterscan allow for the nonzero z p to be small if there is no mRNA demand.Figure 7 shows an example of an ORFM with 10 ribosomes, where the circuit pool z c has a feedbackcontroller. The golden two-codon mRNA is the distinguished x F that produces the inhibitor protein F with λ F = [0 . , , . The comparatively low initiation rate λ F allows other mRNA if present to takeribosomes from z C first. The inhibitor parameters are F p = 4 , γ p = 2 , and δ p = 0 . . Simulations arerun for t = 0 . . . .With no inhibitor and no circuit mRNA, the system in Figure 3a has a total of N Hr = 7 . ribosomesin the host pool and mRNA. With the inhibitor in Figure 7a, the total number of host ribosomes risesto N Hr = 7 . . Figure 3b has a high circuit demand, and N Hr = 3 . . With an inhibitor in Figure 7b, N Hr = 4 . . At no demand, the steady state F p = 0 . , and at high demand, F p = 0 . .V. O PTIMIZATION The goal of optimization in this section is to maximize the rate of protein production by mRNA.Optimization of protein throughput in a single strand has previously been treated by Poker et. al. in [7].Given an n -codon mRNA with allowable choices of λ j in a convex set, finding a λ ∗ that maximizesthe throughput y = λ n e n is a concave optimization problem.Ribosomes with a pool have an additional layer of feedback. For a single strand with λ j taking ribosomesfrom pool z , the effective intake rate as in Equation (1) is λ G ( z ) . Since G ( z ) is monotonically increasingand concave, the throughput y is also a concave function of z .This concavity breaks when multiple strands of mRNA are connected to the same pool, as shown inFigure 8 for an RFMNP with 6 strands (Section II-B). The left plot shows the protein production y i rate of each strand for N r ∈ [0 , . All curves are monotonically increasing with inflection points.This nonconcavity comes from competition effects: an exceptionally greedy strand of mRNA will takeribosomes until it is saturated, and then other strands of mRNA can translate. In contrast, the right handplot illustrates the translation rate as a concave function of tank capacity. The mapping N r → z carriesthe competition effects and nonconcavity. It is conjectured that the mapping N r → y is concave if thereis only a single strand connected to the pool. N r y p Throughput vs. Total Ribosomes z y p Throughput vs. Tank Ribosomes Fig. 8: Protein throughput per mRNAThe objective in this section is to maximize the weighted sum of protein production y w = (cid:80) p,i w pi y pi .Appropriate selection of the weights w can specify particular proteins as desirable or repellent. As aproblem setting, consider a multi-species RFM network with M species and N ribosomes. There willbe M + 1 pools: z E for the empty ribosomes and z p for the translating species. Assume that p = 1 represents the host ribosome with host genes, and that each circuit gene may only accept ribosomesfrom a single pool z p . The matching of pools and circuit genes that maximizes y w requires solvinga binary assignment problem. If a new synthetic gene was to be designed in an environment wherethere is severe competition for circuit ribosomes, it may be favorable for this new gene to employhost ribosomes. This can be generalized to multi-species orthogonal ribosomal networks. An alternativeproblem can be formulated in terms of equilbrium constants K p . If there exists sufficient flexibility toadjust k + p ∈ [ k + ,minp , k + ,maxp ] > and k − p ∈ [ k − ,minp , k − ,maxp ] > , then: K p ∈ [ K minp , K maxp ] = (cid:20) k + ,minp k − ,maxp , k + ,maxp k − ,minp (cid:21) > . For each point K = { K p } Mp =1 , the weighted protein output y w ( K ) can be obtained by finding the steadystate with respect to K using Algorithm 2 and then evaluating y w = (cid:80) w pi y pi . An optimization problemto maximize y w at steady state can therefore be formulated: y ∗ w = max K p (cid:88) p,i w pi y pi , (10) subject to : K p ∈ [ K minp , K maxp ] . (Dynamics in (4a)-(4e).)The search variables are K , and the steady states ( z E , z p , e pij ) are derived quantities of K . Genericsolvers such as a Grid Search, Bayesian Optimization, and Basin Hopping can be used to approximate K ∗ . It can be verified that ∂e pij /∂K p > and ∂e p (cid:48) ij /∂K p < for p (cid:48) (cid:54) = p by evaluating derivatives ofsteady states. Raising K p increases z p , e pij , y pi at the expense of z E and z p (cid:48) , e p (cid:48) ij , y p (cid:48) i for p (cid:48) (cid:54) = p .Fig 9 shows an ORFM with s = 6 and s = 3 with an objective of maximizing total protein output ( w = ) . If all mRNA were connected to a common pool of ribosomes (RFMNP), the resultant output is y w = 3 . . When K = [5 , , then y = 3 . as shown in Figure 9a. Solving the optimization problemin Equation (10) with K min = 1 / and K max = 5 results in K ∗ = [5 , . and the protein output is y w = 3 . in Figure 9b. At steady state in the optimal allocation, z H will contain of all ribosomesin the pools ( ¯ N r ) . The optimization landscape of y w vs. (log ( K ) , log ( K )) is displayed in Figure10. a) Parity: y W = 3 . (b) Optimal: y W = 3 . Fig. 9: Optimize the total protein output ( w = ) -0.6 -0.4 -0.2 0 0.2 0.4 0.6 Log (K ) -0.6-0.4-0.200.20.40.6 Log ( K ) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Fig. 10: Contours of total protein outputVI. C ONCLUSION Competition for finite resources are inevitable in protein translation. Orthogonal ribosomes have beendeveloped to boost protein throughput by decoupling circuit genes from the host pool of ribosomes.We extended the existing RFM to orthogonal ribosomes, and generalized the system to an arbitrarynumber of ribosomal species. Stability results through RLFs and a simple algorithm to compute steadystates were shown for the ORFM system. A self inhibiting feedback controller can adjust ribosomalproduction as needed. Maximizing the weighted sum of protein throughput can be formulated as a low-variable nonconvex optimization problem. Future work includes matching results with lab experimentsand allowing for cross-talk in translation. PPENDIX : P ROOF OF T HEOREM s p = s, n p,i = n , for all p, i . Also, assume m p,i = 1 , for all p, i . The argumentfor the general case will be the similar. Denote σ ( x ) := sgn ( x ) .To use the techniques in [16], [13], we lift (4a)-(4e) to a higher-dimensional space by defining the vacancy w p,ij := 1 − x pij , for all j, i, p . Hence, terms of the form (1 − x pij − ) x pij take the familiar Mass-Action form: w pij − x pij . We will generalize this further by considering arbitrary monotone functions ofthe form R ( w pij − , x pij ) . Hence, we write (4a)-(4e) as: ˙ z E = (cid:80) p R − p ( z p ) − (cid:80) p R + p ( z E )˙ z p = R + p ( z E ) − R − p ( z p ) + (cid:80) i (cid:0) R pin ( x pin ) − R pi ( w pi , z p ) (cid:1) ˙ x pi = R pi ( w pi , z p ) − R pi ( w pi , x pi )˙ x pij = R pij − ( w pij , x pij − ) − R pij ( w pij +1 , x pij ) (11) ˙ x pin = R pin − ( w pin , x pin − ) − R pin ( x pin ) . We only assume that the rates R p,ij , R ± p are monotone w.r.t their reactants, see [13] for the full assump-tions.Consider V ( x, z ) = (cid:80) i,p,j | ˙ x pij | + (cid:80) p | ˙ z p | + | ˙ z E | . We show that V is non-increasing. Using (11), notethat the space can be partitioned into regions, and V is linear in rates on each region. Fix an open region W and write V = (cid:80) p,i,j α pij R pij ( x, z ) + (cid:80) p β p ( R − p ( z p ) − R + p ( z E )) on W . Note V is differentiableon W , and the signs of the currents ˙ x pij , ˙ z p , ˙ z E are constant on W . We claim that the derivative ofeach term in V is non-positive. We show first that α pij ˙ R pij ≤ for all p, i, j . We study three cases: j = 0 , n , < j < n . If j (cid:54) = 0 , n , then R pij appears in ˙ x pij , ˙ x pij +1 . If both have the same sign on W , then(11) implies α pij = 0 . Therefore, α pij (cid:54) = 0 implies σ ( α pij ) = − σ ( ˙ x pij ) = σ ( ˙ x pij +1 ) = − σ ( ˙ w pij +1 ) . Hence, α pij ˙ R pij = α pij ( ∂R pij ∂x pij ˙ x pij + ∂R pij ∂w pij +1 ˙ w pij +1 ) ≤ as claimed, where non-negativity of the partial derivativesfollows from monotonicity. Next, consider the case j = 0 where R pi appears in two currents ˙ x pi , ˙ z p .Similar to the previous case, we conclude that if α pi (cid:54) = 0 then σ ( α pi ) = − σ ( z p ) = − σ ( ˙ w pi ) . Since R pi is a monotone function of w pi , z pi , then α pi ˙ R pi ≤ as claimed. Finally, if j = n , similar analysis can berepeated to conclude that if α pin (cid:54) = 0 then σ ( α n ) = − σ ( x pin ) and α n ˙ R pin ≤ as claimed. Next, we wantto show that β p ˙ R − p ≤ . Fix p . Note that R − p appears in two currents: ˙ z p , ˙ z E . If they have the same signthen β p = 0 . If they have different signs then (11) implies that σ ( β p ) = − σ ( ˙ z p ) = σ ( ˙ z E ) . Since R − p isa function of z p , then β p ˙ R − p ≤ . A similar argument can be replicated to show that − β p ˙ R + p ≤ . Since W and p, i, j were arbitrary, we conclude that ˙ V ≤ over the interior of all regions. The boundariesbetween the regions can be dealt with via the definition of the upper Dini’s derivative, see [16]. Thisconcludes the proof of the first statement of Theorem 1.We proceed to prove part 2. The RLF utilized is non-strict, hence it cannot yield Global AsymptoticStability (GAS) directly. The LaSalle argument proposed in [16] can be used, but it is long. Alternatively,we will show that the Jacobian of (11) (reduced to the invariant space for a fixed N r ) is robustly non-degenerate. This, coupled with the RLF, implies the GAS of any steady-state [17], [13].To prove non-degeneracy, it has been shown in [18], [13] that the Jacobian J of (11) at any pointin the space can be written as J = (cid:80) L(cid:96) =1 ρ (cid:96) Q (cid:96) for some ρ (cid:96) ≥ , rank-one matrices Q (cid:96) and some L . Here, L = M (2 + s (2 n + 1)) . The coefficients { ρ (cid:96) } correspond to the partial derivatives of therates { R pij } , { R ± p } , while { Q (cid:96) } correspond to the network structure and are independent of the rates.Furthermore, it has been shown also [13] that the existence of an RLF implies that it sufficient to find ne positive point ρ ∗ = ( ρ ∗ , .., ρ ∗ L ) for which the (reduced) Jacobian is non-degenerate to show that isnon-degenerate for all ρ (cid:96) > . We will find that point next by studying the structure of the Jacobian.Similar to a single pool [11, SI], it can be easily seen that J has non-negative off-diagonals and strictlynegative diagonals (i.e, J is Metzler). In addition, conservation of the number of ribosomes impliesthat T J = 0 . By the definition of the Jacobian, all the entries in each column contain only thepartial derivatives with respect to the state variable associated to the column. Hence, we can choose thecorresponding ρ (cid:96) ’s such the diagonal entry in each column is scaled to − . Therefore, we consider theJacobian evaluated at the chosen point ρ ∗ such that J ∗ = P − I , where I is the identity matrix and P a nonnegative irreducible column-stochastic matrix. By Perron-Frobenius Theorem, P has a maximaleigenvalue 1 with algebraic multiplicity 1. Therefore, J ∗ has a single eigenvalue at and the remainingeigenvalues have strictly negative real-parts. Hence, the reduced Jacobian at ρ ∗ is non-degenerate. Robustnon-degeneracy and GAS follows.The existence of a steady-state follows from Brouwer’s fixed point theorem since (11) evolves in acompact space (for a fixed N r > ), and uniqueness follows from non-degeneracy and GAS. Thepositivity of the steady-state follows from persistence of the ORFM which can be shown graphicallyby the absence of critical siphons [19]. (cid:4) Acknowledgements: This research was partially funded by NSF grants 1849588 and 1716623. JaredMiller thanks Prof. Mario Sznaier for his support and discussions.R EFERENCES [1] O. Rackham and J.W. Chin. A network of orthogonal ribosome · mRNA pairs. Nature Chemical Biology , 1(3):159–166, 2005.[2] A. Costello and A.H. Badran. Synthetic biological circuits within an orthogonal central dogma. Trends in Biotechnology , In press,2020.[3] L.M. Chubiz and C.V. Rao. Computational design of orthogonal ribosomes. Nucleic acids research , 36(12):4038–4046, 2008.[4] C.T. MacDonald, J.H. Gibbs, and A.C. Pipkin. Kinetics of biopolymerization on nucleic acid templates. Biopolymers: OriginalResearch on Biomolecules , 6(1):1–25, 1968.[5] C.T. MacDonald and J.H. Gibbs. Concerning the kinetics of polypeptide synthesis on polyribosomes. Biopolymers: Original Researchon Biomolecules , 7(5):707–725, 1969.[6] S. Reuveni, I. Meilijson, M. Kupiec, E. Ruppin, and T. Tuller. Genome-scale analysis of translation elongation with a ribosomeflow model. PLoS Comput Biol , 7(9):e1002127, 2011.[7] G. Poker, Y. Zarai, M. Margaliot, and T. Tuller. Maximizing protein translation rate in the non-homogeneous ribosome flow model:a convex optimization approach. Journal of The Royal Society Interface , 11(100):20140713, 2014.[8] Y. Zarai, M. Margaliot, E.D. Sontag, and T. Tuller. Controllability analysis and control synthesis for the ribosome flow model. IEEE/ACM Trans. Comput. Biol. Bioinf. , 15(4):1351–1364, 2017.[9] M. Margaliot and T. Tuller. Ribosome flow model with positive feedback. J. R. Soc., Interface , 10(85):20130267, 2013.[10] A. Raveh, Y. Zarai, M. Margaliot, and Tamir Tuller. Ribosome flow model on a ring. IEEE/ACM transactions on computationalbiology and bioinformatics , 12(6):1429–1439, 2015.[11] A. Raveh, M. Margaliot, E.D. Sontag, and T. Tuller. A model for competition for ribosomes in the cell. Journal of The RoyalSociety Interface , 13(116):20151062, 2016.[12] A.P.S. Darlington, J. Kim, J.I. Jim´enez, and D.G. Bates. Dynamic allocation of orthogonal ribosomes facilitates uncoupling ofco-expressed genes. Nature Communications , 9(1):695, 2018.[13] M. Ali Al-Radhawi, D. Angeli, and E.D. Sontag. A computational framework for a Lyapunov-enabled analysis of biochemicalreaction networks. Plos Computational Biology , 16(2):e1007681, 2020.[14] M. Ali Al-Radhawi and D. Angeli. Piecewise linear in rates Lyapunov functions for complex reaction networks. In Proceedings ofthe 52nd IEEE Control and Decision Conference (CDC) , pages 4595–4600, 2013.[15] F. Blanchini and G. Giordano. Piecewise-linear Lyapunov functions for structural stability of biochemical networks. Automatica ,50(10):2482 – 2493, 2014.[16] M. Ali Al-Radhawi and D. Angeli. New approach to the stability of chemical reaction networks: Piecewise linear in rates Lyapunovfunctions. IEEE Trans. on Automatic Control , 61(1):76–89, 2016.[17] F. Blanchini and G. Giordano. Polyhedral Lyapunov functions structurally ensure global asymptotic stability of dynamical networksiff the Jacobian is non-singular. Automatica , 86:183–191, 2017.[18] M. Ali Al-Radhawi and D. Angeli. Robust Lyapunov functions for complex reaction networks: An uncertain system framework. In Proceedings of the IEEE 53rd Conference on Decision and Control (CDC) , pages 3101–3106, Dec 2014.[19] D. Angeli, P. De Leenheer, and E.D. Sontag. A Petri net approach to the study of persistence in chemical reaction networks.