A. Background IVING cells have the ability to perform sophisticated operations that include maintaining homeostasis against noise, responding appropriately to various input signals, constructing complex structures such as proteins, and adapting to novel environments. Reverse engineering the biochemical circuits responsible for implementing such operations has revealed various control mechanisms that include regulation of gene expression via transcription factors (TFs) and/or noncoding RNAs [1], [2]. This has inspired the development of engineering approaches that mimic natural circuits by inserting new synthetic circuits into cells to modify their behavior or create new functionalities. Applications are wide-ranging and include immunotherapy, programmed micro-organisms for diagnostics and therapy, biofuel production, and many others [3], [4], [5]. Despite the great promise, numerous challenges exist. In particular, genetic circuits utilize common resources for transcription and translation such as RNA polymerases (RNAPs), ribosomes, tRNAs, and others. Insertion of new circuits increases the load on the cell's reservoirs. This, in turn, can create indirect interactions [6] that impede the proper functioning of the circuit, retard cellular growth, or lead to premature apoptosis [7], [8]. Several approaches have been proposed to ameliorate this problem, including dynamic control [9], orthogonal ribosomes and RNAPs [10], [11], [12], forward engineering of the circuit to account for resource competition [13], and distributed computation [14], [15].
Several of the aforementioned approaches assume that it is possible to identify the mode of competition and the limited resources responsible for performance degradation. However, it is not always possible to infer the correct model of competition from the expression data of the circuit. Possible competition effects to account for include promoters competing for RNAPs, mRNAs competing for ribosomes, transcription factors competing for promoters, enzymes competing for substrates, substrates competing for enzymes, etc.
In this work, we identify competition phenotypes, i.e., features that could allow one to distinguish the scarce resource (or resources) responsible for performance deterioration. Are there qualitatively different types of competition effects? Are there equivalent effects that can be treated in a unified manner? Answers to these questions will help guide theoretical analysis as well as the design of targeted interventions that mitigate undesirable effects through, for example, the use of feedback control to regulate the level of a scarce resource, or the optimization of appropriate circuit parameters. We will ask how different aspects of gene expression are impacted by two factors: (1) the level of the resource being sharedsuch as an activating or repressing transcription factor (TF), RNAPs, or ribosomes-and (2) the level of competition from other biochemical species -such as other genes or mRNAs. We describe interventions that can be used to distinguish the different competition phenotypes. We also discover instances where, conversely, competition for different resources might result in mathematically equivalent competition prototypes.
Notation: Chemical species are denoted by non-italic large caps, time-dependent state variables are denoted by small caps, and constants are denoted by large cap italics. For example, a protein is denoted by Y , its time-dependent concentration is denoted by y(t) while its steady-state value is denoted by Y .
1) The system: Consider a synthetic circuit with an external input U and internal state vector x which includes the concentrations of promoters, mRNAs, proteins, etc. The output 0 P n i 0
Fig. 1.
Framework for studying the problem of competition model discrimination.
is denoted by y. The circuit utilizes the free limited resource r which is also utilized by other competing (or interfering) circuits with inputs I1 ; ::; In and internal state vectors z1; ::; zn. Figure 1 provides a pictorial representation of the system.
We model the entire system by a system of ordinary differential equations structured as follows:
x _ = f (x; r; U );
(1) z _i = gi (zi ; r; Ii ); i = 1; ::; n;
(2) r_= '(r; x; z1 ; ::; zn );
for some C 1 vector fields f ; g1; ::; gn ; ', and function h. We assume that the limited resource is conserved. It partakes in the synthetic and competing circuits without being consumed or annihilated. Concretely, there exists nonnegative vectors d0; ::dn of compatible dimensions such that r + d T x + i = 1 d T zi = R T ;
(5) where R T is the total resource, and r is the free (available) resource. In order for (1),( 2),(3) to satisfy (5), the following relation is assumed to be satisfied: r_ + d T x _ + i = 1 d T z _i 0. In this paper, we perform our analysis at steady-state. We will assume that for each choice of the inputs U; I1; ::; In and total resource R T , there exists a steady-state (X ; Z1 ; ::; Zn ; R) that is globally asymptotically stable.
After eliminating the intermediate variables, the corresponding steady-state output Y can be written as: Y = H ( R T ; U; I1; ::; In); (6) for some function H .
2) Performance evaluation: The performance of a circuit with competition is compared to its performance with no competition. For this purpose, the Competition-induced Performance Deterioration Ratio (CPDR) is defined as: : = Y = Y j competition=0 ; where the total competition is I : = i = 1 Ii . 3) Problem Formulation: The formulation in (1)-( 4) implicitly assumes that r is the limited resource, while other resources are abundant and appear as kinetic constants in the functions f; g1; ::; gn; h. Changing the identity of the limited resource will change the model (1)-(4). Our aim is to compare the qualitative differences in the steady-state input-output data that follow from the scarcity of different resources.
The paper is organized as follows. In section II, we discuss transcription/translation systems where the mRNA is protected by the ribosomes, while we drop that assumption in section III. We conclude with a discussion in section IV. The input to the circuit is the total concentration of the input promoter U which is assumed to be constant; the unbound (free) promoter is denoted by Uf .
We first assume that mRNAs are protected from decay when bound to ribosomes, which is a valid assumption in many situations [16]. The transcription and translation reactions can be then written as follows:
where Q denotes RNAP, R denotes the ribosome, M denotes mRNA, E is the promoter-RNAP complex, and F is the mRNA-ribosome complex. Free competing promoters If ; i ; i = 1; ::; n can bind to RNAP and produce mRNAs Mi which can bind to the ribosome. This can be written as:
The promoters are conserved. Hence, we have:
I f ; i + E i = I i ; i = 1; ::; n:
Let K j = k j =kj ; Aj = j = j ; Wj = w j =wj ; B = = , j = 1; 2 be RNAP Dissociation Ratio (DR), tran-scription ratio, ribosome DR, and protein expression ra-tio, respectively.
By writing the irreversible reactions as ;
M, ; ) Y ; ; ) M i , ;
F i Y i , the steadystates of 1 the network ( 7)-( 8),( 9)-( 10) can be computed by noting that the network is detailed-balanced, i.e. the forward and backward rates in each reaction are equal at steady state. Hence, the steady state values of the mRNAs M; Mi are determined by the equilibrium values A 1 E ; A 2 E i ; i = 1; ::; n, respectively, and are independent of the translation process. We get the following steady-state expressions:
Note that Q; R are the levels of the free RNAP and ribosome, respectively.
In the formalism defined in §1.B and depending on the identity of the limited resource, Eq. ( 1) describes the dynamics of ( 7)-( 8), while Eq. ( 2) describes the dynamics of ( 9)- (10). In order to write the output in the form ( 6), we will study several scenarios in which either the ribosome is limited, or RNAP is limited. We keep the other resource abundant to isolate the effects of the limitations in a single resource.
Remark 1: A slightly different mathematical model of translation (8) would have the mRNA/ribosome complex dis-
This article has been accepted for publication in IEEE Control Systems Letters. This is the author's version which has not been fully edited and Remark 2: An alternative way to regulate a target circuit is to use a small molecule such as A H L to activate the target promoter [7]. The RNAP binding reactions can be written as follows (compare to ( 7)): U ) Uf ; Uf + Q * E, where U is the free inactive promoter, and is proportional to the exogenous input. Our model encompasses this case also since it can be shown that the reparameterization K 1 = (k 1 =k1 )(1 + =) recovers the equations ( 13)- (15). In other words, the effective effect of the external input in our model is to regulate the RNAP dissociation ratio K 1 . Note that A H L concentrations can be precisely regulated over four orders of magnitude [8]. Therefore, employing promoters inducible by A H L will not have a considerable effect on the dynamic range of K 1 .
Parameter Ranges: In order to keep the numbers within biological ranges, we use the following parameters: K i 2 [0:3; 10000] nM; Wi 2 [5; 2000] [7], [8]. Furthermore, whenever the following parameters are not treated as variables, we use the following numbers as in [7]: QT = 500 nM; RT = 1000 nM; Ai = 10; B = 300, i = 1; 2.
In this case, RNAP Q is unaffected by the circuit ( 7)-( 8),( 9)-( 10), hence its level Q will be constant, Q = QT .
Let F I = i = 1 Fi , then the conservation law (5
Solving the resulting algebraic equations, the free resource at steady state is found to be:
The output is linearly dependent on R T , and takes an inhibiting Michaelis-Menten form with respect to the competition I . It also takes a Michaelis-Menten form with respect to U . The CPDR is given as: = (1 + 1 U )=(1 + 1 U + 2 I ); which is independent of the total resource R T .
We next study the properties of the output as it depends on the input and the competition.
As a function of the resource R T , the output is linear as noted above, but is concave with respect to U . As a function of the competition, it is convex. Note d 2 Y =dI 2 > 0 for all I > 0 . Figure 2 depicts the typical behavior of the competition phenotype associated with this mode of competition.
Remark 3: The fact that the CPDR is independent of the total resource might lead one to think that the total resource is irrelevant for reducing competition. However, this depends on how we define competition reduction. If we consider strategies to reinstate the level of the output Y to its competition-free level Y j , then we can increase the total resource to compensate for the reduction in the output due to competition. In particular, we can write the required total resource as R = 1 + 1 U + 2 I Y j I = 0 .
In this case, the ribosome R will be unaffected by the circuit ( 7)-( 8),( 9)-( 10), hence its level R will be assumed to be constant. In other words, we have R = R T . Therefore, we write the conservation law (5) as:
In the most general case, solving (17) for Q requires solving a cubic equation. Therefore, in this subsection, we assume that K 1 = K 2 = K to simplify the analysis. This can be justified in cases when the gene and its competitors are located on the same plasmid. For instance, the RNAP DRs are taken to be K 1 = K 2 = 200 nM when simulating the data in [7]. In this case, the free RNAP Q and the output Y are
As QT grows without bound, we have limQ ! 1 Y = 1 U R T , i.e., the protein is expressed at maximum capacity. If I grows without bound, then l i m I ! 1 Y = 0.
The competition reduction ratio is written as follows:
which depends on the total resource unlike the case in §2.1. Convexity of the output as a function of the resource: The output is globally concave, because:
for all QT 0. The typical competition phenotypes are plotted in Figure 3-a,b. Convexity of the output as a function of the competition : Simulations show that the response is globally convex when QT is small. For larger QT , the response starts concave and then it has an inflection point. Using the same parameters above, when QT = 1 the response is convex at zero as verified by computing the second derivative of Y with respect to I at zero. Figure 3-b shows the transition from convexity to concavity with higher QT . For both panels we have
C. Distinguishing the competition phenotypes a) Convexity/concavity: We have been able to prove that, at steady state, limited ribosome and limited RNAP result in qualitatively distinct competition phenotypes. For the second, the output takes a Michaelis-Menten form and is globally convex (18) with respect to the resource (RNAP) level as shown in Figure 3-a, but, in contrast, for the first it is perfectly linear with respect to the resource (ribosome level) regardless if it is abundant or limited; see Fig. 2-a. Furthermore, a high level of RNAP provides buffering against low levels of competition in the second case as shown in Figure 3-b, while the output drops quickly even with low competition in the first case; see Fig. 2b. We are able to characterize this mathematically by proving that the output is initially concave (i.e., superlinear) with respect to competition mediated by limited RNAP, while it is globally convex (i.e., sublinear) with limited ribosomes. Thus, using either criterion (titrating resource, or titrating competitors), we can see a clear difference between these two types of context limitations. For instance, the level of competition can be controlled by adjusting the dosage of an inserted plasmid [7]. Hence, our criteria can be checked visually from the plot of the output versus the competitor dosage. From a different point of view, theory helps us guess the source of competition based upon experimental data.
b) Effect of the RNAP dissociation ratio: The effective RNAP DR K 1 can be modified by using an inducible promoter (see Remark 2). Keeping the remaining parameters fixed, we may think of the outputs Y and YI : = n Yi as functions of K 1 . Let us now analyze the trade-off between the two outputs when varying the parameter K 1 while keeping the resource levels constant. The resulting parametric curves are also known as isocost curves in the language of economics [7]. Such relations can be derived by solving Eq. ( 16) for K 1 in terms of Y , and then substituting it in the expression of YI .
In the case of L R A P (i.e., Q = QT ), it can be shown that the relationship between the parameterized outputs is linear and is given by the following parametric equation:
Note that ( 19) is independent of the ribosome DR W1.
In comparison, the case of L PA R (i.e., R = R T ) is more complicated. The corresponding relationship between Y I ( K 1 ) and Y ( K 1 ) is generally nonlinear and its computation requires solving a cubic equation, as noted before. In order to probe the effect of K 1 , we use the simplifying assumption Q K 1 ; K 2 which is often satisfied by practical systems [7]. Under this approximation, it can be shown again that the relationship between Y ( K 1 ) and Y I ( K 1 ) is linear and is given by the parametric equation:
In this case, the relationship depends on W1. Therefore, the two modes of competition (for RNAP or ribosomes) are in principle distinguishable by modifying W1 via variable Ribosome Binding Site (RBS) strengths. Note that this conclusion stands even without the assumption made earlier.
Even though the linear approximation (20) may not hold in all situations, dependence on W1 in the K1-parameterized relationship between Y and YI indicates limited RNAP as seen in Figure 4-a) which depicts a simulation of the parameterized curves for various values of W1. In comparison, Eq. ( 19) is linear and is independent of W1. Figure 4-b) shows experimental data from [7] that depicts the same scenario. Our theoretical and computational prediction is consistent with the slope change noticed in the experimental data. We derived our conclusions assuming that RNAP is limited and that ribosomes are abundant. In practical situations, both resources might be in limited supplies [7].
c) Effect of the Ribosome Dissociation Ratio: Let us now consider modifying W1, instead of K 1 . We can again write the outputs as Y (W1); YI (W1). In the case of LRAP, the relationship is again linear and it can be written in the same form as (19). It can be seen that it is independent of the RNAP DR K . In the case of LPAR, under the assumption Q K 1 ; K 2 , we get: YI (W1 ) = A 2 B I K 1 Q T R T ; which depends on the RNAP DR K 1 but is constant. This conclusion holds without the assumption above since the free RNAP Q which solves (17) is independent of ribosome DR W1. Hence, the expressions Y (W1); YI (W1) are not related.
Figure 5-a) shows that the case of limited ribosomes manifests as a decreasing linear relationship when parameterized by W1 (which can be experimentally controlled by R B S strengths), but it is independent of K 1 (which can be controlled by utilizing an inducible promoter). On the other hand, the same relationship is constant in the case of limited RNAP, but the level is dependent on K 1 . This manifests by examining the YI -axis intercepts, which depend on K 1 in the case of LPAR, but are independent of it in the case of LRAP. d) Effect of the total copy number of the competitor: We consider next modifying I while keeping U fixed. We write the outputs as Y (I ); YI (I ).
In the case of LRAP, we get the linear relationship: Y I ( I ) = R T W 1 ( K 1 + Q T ) + 1 Y (I ); where the YI -axis intercept depends only on R T .
In the case of LPAR, the relationship is generally nonlinear, but under the assumption Q K 1 ; K 2 we get the following linear relationship: Y
Y (I ); where the YI -axis intercept depends on the ribosome DR W1. Figure 5-b) depicts the parametric curves corresponding to L R A P case where it can be seen that the YI -axis intercept is constant for various W1.
e) Summary: The results are summarized in Table I. [16]. To model this, we modify ( 7)-( 8),( 9)-( 10) by adding the following reactions which describe the decay of mRNA while bound to the ribosome: F ! R ; F i ! R ; i = 1; ::; n:
(21
The network is no longer detailed balanced after the inclusion of decay of the bound mRNA-ribosome complexes. The steady-values of M , Mi are now dependent on the resource R . But we still have the promoter-RNAP complex expressed at steady-state as:
Let aj = = j be the ratio of the decay rates of the bound and unbound mRNA. The steady-state values of the mRNAs are given as: M = A1 E =(1 + Ra1=W1); Mi = A2 Ei =(1 + Ra2=W2); where Wj = ( + w j )=wj ; j = 1; 2. The steady-values of the mRNA-ribosome complexes and the output in terms of the free resource R are given as:
Observe that when = = 0; we recover the case discussed in §2. We study different scenarios next.
As before, we have abundant RNAP, hence Q = QT . The only conservation law is (16). Therefore, we need to solve for the free ribosome R . Let E be as given in ( 22), and let E I = E i where E i is defined in (22). The conservation law (16) results in a cubic equation. Therefore, we assume that A = A1 = A2; W = W1 = W2 ; K = K 1 = K 2 ; a = a1 = a2 in order to simplify the analysis. We get (23). The output can be written as in (24). Note that when a1 = a2 = 0, we get the case discussed in §2.1. Next, we show that the properties of the system above can be deduced by studying a different system that has been studied earlier.
It is perhaps surprising, that two very different biological systems may lead to mathematically identical competition phenotypes. To formalize this, let us write the (6) as: Y = H (V ; ) where the inputs are V = [U; RT ; I1 ; ::; In ] T , and are a set of parameters (kinetic rates, for example) that appear in (1)-(4). Hence, we have a steady-state expression H1(V1; 1 ) that gives us the amount of output, as well as a second function H2(V2; 2 ) for a different system; equivalent phenotypes will have the property that there is a diffeomorphism : (V1; 1 ) ! (V2; 2 ) so that every function H2(V2; 2 ) can be written as H1 ( (V1; 1 )).
As a concrete example, let us revisit scenarios §2.2 and §3.1 discussed earlier with A1 = A2 = A. For the system in §3.1, we consider the case in which the mRNA decays at the same rate, whether it is bound to the ribosome or not, i.e., a1 = a2 = 1. One can prove that the two systems are equivalent, under a reparameterization given as follows: k =k ! (w + )=w; QT ! R T ; U ! U QT =( ( K + QT )); B A U R T =(W + R T ) ! B ; I ! I Q T =( ( K + QT )); where I : = Ii . Note that the underlying biochemical systems are very different, and the two systems of ODEs are distinct. In fact, they result in different transient behavior. However, the steady-states are the same, as shown theoretically and illustrated in Fig. 6.
Similar to the previous section, we let R = R T . Hence, the only conservation law is (17). Since decay affects only the with the rest of the genome [10], [7], [11]. Needless to say, this does not always hold. In addition, we have assumed a fixed number of competitors and a constant amount of allocated resources. However, the number of active pathways and the amount of allocated resources in a cell change dynamically depending on stress and growth conditions [19]. Studying such scenarios is subject to future work. ribosome-mRNA complex, we can see immediately that this case is very similar to the case discussed in §2.2 except for an additional factor. Hence, similar analysis can be replicated.
1) Generalization: The competition phenotypes studied in the previous sections can be generalized to other biologi-cal contexts by classifying them into two main categories: externally-regulated targets, and conserved targets. Ribosome competition studied in §2.1 can be studied under the first category. This is since the target which consumes the resource is the mRNA which is not conserved. Hence, the model can be essentially written as follows: In addition to mRNAs, the above model can represent ligandactivated enzymes M; Mi; i = 1; ::; n competing for a substrate R, single-guide RNAs (sgRNAs) competing for a limited amount of dCas9 in CRISPRi [17], [18], or externallyactivated TFs competing for a single promoter. On the other hand, the case discussed in §2.2 has the RNAP as a limited resource. Hence, the target that consumes the resource is the promoter which is conserved. Hence, the model can be written essentially as follows: In addition to promoter, the above model can represent conserved enzymes M; Mi; i = 1; ::; n competing for a substrate R, or conserved TFs competing for a single promoter.
2) Limitations: Our framework has included multiple simplifications to allow analytical derivations and facilitate clear interpretations. For instance, we assumed that the competitors behave similarly to each other. This can be justified when the target gene and its competitors are co-located on the same plasmid with high copy numbers, and with their own orthogonal RNAPs and ribosomes to minimize interactions
© 2022 IEEE. Personal use is permitted, but republication/redistribution requires IEEE permission.See https://www.ieee.org/publications/rights/index.html for more information.