Integral feedback in synthetic biology: Negative-equilibrium catastrophe
IIntegral feedback in synthetic biology:Negative-equilibrium catastrophe
Tomislav Plesa, ∗ Guy-Bart Stan, Thomas E. Ouldridge, Alex Dack Abstract : A central goal of synthetic biology is the design of molecular controllers that can manip-ulate the dynamics of intracellular networks in a stable and accurate manner. To address the factthat detailed knowledge about intracellular networks is unavailable, integral-feedback controllers(IFCs) have been put forward for controlling molecular abundances. These controllers can maintainaccuracy in spite of the uncertainties in the controlled networks. However, this desirable featureis achieved only if stability is also maintained. In this paper, we show that molecular IFCs sufferfrom a hazardous instability called negative-equilibrium catastrophe (NEC), whereby all nonneg-ative equilibria vanish under the action of the controllers, and some of the underlying molecularabundances blow up. We show that unimolecular IFCs do not exist due to a NEC. We then de-rive a family of bimolecular IFCs that are safeguarded against NECs when unimolecular networksare controlled. However, we show that, without a detailed knowledge of the controlled systems,NECs cannot be prevented when IFCs are applied on bimolecular, and hence the vast majority ofintracellular, networks.
A main objective in synthetic biology is achieving control over living cells [1, 2, 3, 4, 5, 6]. Thischallenging problem requires manipulation of the intricate and only partially known intracellularmachinery, involving nonlinear and noisy biochemical interactions [7, 8, 9]. In this context, the taskis to design a molecular controller network that, when embedded into a black-box input network,guarantees that the resulting output network autonomously traces a predefined dynamics in a stable and accurate manner. Control can be sought over deterministic dynamics when all of the molecularspecies are in higher-abundance [10], or over stochastic dynamics when lower copy-number speciesare present [31]. See also Figure 1, and Appendices A and B.Molecular abundances and parameter values of intracellular networks change over time, andcan be experimentally measured only with limited precision. To overcome these challenges, con-trollers must ensure that both stability and accuracy are maintained robustly with respect to theuncertainties in the initial conditions and parameter values of the input networks. When successful,such controllers ensure robust adaptation (homeostasis) - a fundamental design principle of livingsystems [12, 13, 14, 15, 16, 17]. Furthermore, it is of critical importance that controllers ensurestability and accuracy for generic members in a given class of practically-relevant input networks.A controller that operates as intended for only a negligible subset of input networks is not only oflimited applicability, but also potentially hazardous. In particular, when applied outside its scope,a controller can lead to an unbounded increase of a molecular abundance that can place a lethalburden on the input system [18]. Therefore, it is of great importance to theoretically determinethe set of admissible input networks for which a given controller operates safely.In context of electro-mechanical systems, accuracy can be achieved robustly via so-called integral-feedback controllers (IFCs) [19]. Loosely speaking, given a suitable system, IFCs dynamically cal-culate a time-integral of a difference (error) between the target and actual values of the controlledvariable. The error is then used to decrease (respectively, increase) the controlled variable when itdeviates above (respectively, below) its target value via a negative-feedback loop. However, IFCsthat are implementable with electro-mechanical systems are not necessarily implementable with Department of Bioengineering, Imperial College London, Exhibition Road, London, SW7 2AZ, UK ∗ Corresponding author and lead contact. E-mail: [email protected] a r X i v : . [ q - b i o . M N ] F e b igure 1: Schematic representation of biochemical control . A black-box input network R α is dis-played, consisting of unknown biochemical interactions shown in grey, where X → X (respec-tively, X (cid:97) X ) indicates that species X influences species X positively (respectively, negatively).The particular input network contains one target species X , shown as a green hexagon, that canbe interfaced with a controller. The rest of the input species, that cannot be interfaced with acontroller, are called residual species; one of the residual species, X , is highlighted as a yellowtriangle, while some of the other ones are displayed in grey. A controller network R β,γ = R β ∪ R γ is also shown, consisting of controlling species Y and Y shown as a red circle and blue square,respectively, whose predefined biochemical interactions are shown in black. The controller consistsof the core R β = R β ( Y , Y ) , specifying how the controlling species interact among themselves, andthe interface R γ = R γ ( X , Y , Y ) , specifying how the controlling species Y and Y interact withthe target species X . The composite network R α,β,γ = R α ∪ R β,γ is called an output network. Theparticular controller R β,γ displayed corresponds to the network (8) , with i = j = 1 , from Section 3 . biochemical reactions [20]. Central to this problem is the fact that the error takes both positiveand negative values and, therefore, cannot be directly represented as a nonnegative molecular abun-dance. In this context, a linear non-biochemical IFC has been mapped to a bimolecular one in [21],which has been adapted in [22] and called the antithetic integral-feedback controller (AIFC).Authors from [22] have analyzed performance of the AIFC in context of controlling averagespecies copy-numbers at the stochastic level. In this setting, a number of abstract mathematicalresults have been presented that are valid under a set of technical assumptions ensuring stability.For example, in [22, Theorem 2], the authors specify a class of unimolecular input networks thatcan be successfully controlled with the AIFC; in particular, these input networks have to satisfyan algebraic condition, given as [22, Equation 7]. The authors then put forward a gene-expressionsystem as the input network, given by [22, Equation 9], and demonstrate that the AIFC canarbitrarily control the average protein copy-number. However, this success relies on using a gene-expression input network with no basal level of transcription - an unjustified restriction that isnecessary for the system to always satisfy the conditions of the theorem. Similar choice of theinput network is put forward in [23], where the authors experimentally implement a part of theAIFC to control a protein species with no basal production. The AIFC has also been analyzed incontext of controlling species concentrations at the deterministic level in [24, 25]. However, similarto the running example [22, Equation 9], all of the results in [24, 25] are derived for a severelyrestricted class of unimolecular input networks that always satisfy [22, Theorem 2]. Questionsof critical importance arise: Can the AIFC successfully control practically-relevant unimolecularnetworks? If the control fails, are the consequences biochemically safe or hazardous? Does thereexist a molecular IFC with a better performance than the AIFC? How do IFCs perform whenapplied to practically-relevant bimolecular networks?2a) (c) (b) (d)Figure 2: Caricature representation of a successful and catastrophically failed intracellular control .Panel (a) displays a cell successfully controlled with the IFC in the setup shown in
Figure 1 . Thetime-evolution of the underlying species concentrations are shown in panel (b) . In particular, thetarget species X approaches a desired equilibrium, shown as a black dashed line, and the equilibriumfor the residual species X is positive. Panel (c) displays a cell that has taken lethal damagedue to a failure of the IFC . In particular, as shown in panel (d) , the target equilibrium for X enforces a negative equilibrium for the residual species X . However, since molecular concentrationsare nonnegative, this equilibrium cannot be reached and, therefore, control fails. Furthermore,the failure is catastrophic, as concentrations of some of the underlying species (in this example,species X and Y ) blow up, placing a lethal burden on the cell. Panels (b) , and (d) , are obtainedby solving the reaction-rate equations for the output network (16) ∪ (19) from Section 5 with thedimensionless coefficients ( α , α , α , α ) = (1 , , / , / , ( β , β , γ , γ , γ ) = (100 , , , , ,and with ( α , α , α , α ) = (1 , , / , / , ( β , β , γ , γ , γ ) = (100 , , , / , , respectively. negative-equilibrium catastrophe (NEC),which we outline in Figure 2.The paper is organized as follows. In Section 2, we prove that unimolecular IFCs do notexist due to a NEC. In Section 3, we use the theoretical framework from [20] to derive a classof bimolecular IFCs given by (8). In particular, as a consequence of demanding in the derivationthat the controlling variables are positive, we obtain an IFC that influences the target speciesboth positively and negatively, in contrast to the AIFC that only acts positively. In Section 4,we apply different variants of these controllers on a unimolecular gene-expression network (9), andthen generalize the results to any stable unimolecular input network. We show that the AIFCfrom [22] can lead to a NEC when applied to (9), both deterministically and stochastically. Morebroadly, we prove that the AIFC does not generically operate safely when applied to unimolecularnetworks. We also prove that there exists a variant of the IFCs (8), involving both positive andnegative interfacing, that achieves stability and accuracy for all stable unimolecular input networks.However, in Section 5, using as an example network (19), we demonstrate that the safety profile ofmolecular IFCs cannot be guaranteed when applied to bimolecular input networks. We conclude bypresenting a summary and discussion in Section 6. Notation and background theory are introducedas needed in the paper, and are summarized in Appendices A and B. Rigorous proofs of the resultspresented in Sections 2 and 4 are provided in Appendices C–D and E, respectively. In this section, we consider an arbitrary one-species black-box input network R α = R α ( X ), where X is a single target species and there are no residual species, see Figure 1 for a more general setup.In this paper, we assume all reaction networks are under mass-action kinetics [10] with positivedimensionless rate coefficients, which are displayed above or below the reaction arrows; we denotethe rate coefficients of R α by α ∈ R a> , where R > is the space of positive real numbers. In whatfollows, we say that a reaction network is unimolecular (respectively, bimolecular) if it contains areaction with one (respectively, two), but not more, reactant molecules. Linear non-biochemical controller . Let us consider a controller formally described by thenetwork ¯ R β,γ = ¯ R β ( ¯ Y ) ∪ ¯ R γ ( X , ¯ Y ), given by¯ R β : ∅ β −→ ¯ Y , ¯ R γ : X γ −→ X − ¯ Y , ¯ Y γ −→ ¯ Y + X . (1)Here, ¯ R β = ¯ R β ( ¯ Y ) is the controller core, describing the internal dynamics of the controlling species¯ Y , while ¯ R γ = ¯ R γ ( X , ¯ Y ) is the controller interface, specifying interactions between ¯ Y and thetarget species X from the input network, see also Figure 1. Let us denote abundances of species { X , ¯ Y } from the output network R α ∪ ¯ R β,γ at time t ∈ R ≥ by ( x , ¯ y ) = ( x ( t ) , ¯ y ( t )) ∈ R . At4he deterministic level, formal reaction-rate equations (RREs) [10] readd x d t = f ( x ; α ) + γ ¯ y , x ∗ = β γ , d¯ y d t = β − γ x , ¯ y ∗ = − γ − f (cid:18) β γ ; α (cid:19) , (2)where f ( x ; α ) is an unknown function describing the dynamics of R α , and ( x ∗ , ¯ y ∗ ) ∈ R is theunique equilibrium of the output network, obtained by solving the RREs with zero left-hand sides.Assuming that ( x ∗ , ¯ y ∗ ) is globally stable, network (1) is an IFC; in particular, in this case, x ∗ =( β /γ ) is independent of the initial conditions and the input coefficients α . However, controller (2)cannot be interpreted as a biochemical reaction network. In particular, the term ( − γ x ) in (2)induces a process graphically described by X γ −→ X − ¯ Y in (1), which consumes species ¯ Y evenwhen its abundance is zero. Consequently, ( x , ¯ y ) may take negative values and, therefore, cannotbe interpreted as molecular concentrations [20]. Unimolecular controllers . The only unimolecular analogue of the IFC (1), that containsonly one controlling species Y , is of the form R β : ∅ β −→ Y , R γ : X γ −→ X + Y ,Y γ −→ Y + X . (3)The RREs and the equilibrium for the output network R α ∪ R β,γ are given byd x d t = f ( x ; α ) + γ y , x ∗ = − β γ , d y d t = β + γ x , y ∗ = − γ − f (cid:18) − β γ ; α (cid:19) . (4)Given nonnegative initial conditions, variables ( x , y ) from (4) are confined to the nonnegativequadrant R ≥ , and represent biochemical concentrations. However, the x -component of the equi-librium from (4) is negative and, therefore, not reachable by the controlled system. Furthermore, y is a monotonically increasing function of time, d y / d t >
0, i.e. y blows up. We call this phe-nomenon a deterministic negative-equilibrium catastrophe (NEC), see also Appendix A. Network (3)not only fails to achieve control, but it introduces an unstable species and is, hence, biochemicallyhazardous. In Appendix C, we prove that a NEC occurs at both deterministic and stochastic levelsfor any candidate unimolecular IFC, which we state as the following theorem. Theorem 2.1.
There does not exist a unimolecular integral-feedback controller.Proof.
See Appendix C.
Theorem 2.1 implies that only bimolecular (and higher-molecular) biochemical networks may exertintegral-feedback control. An approach to finding such networks is to map non-biochemical IFCsinto biochemical networks, while preserving the underlying integral-feedback structure. This taskcan be achieved using special mappings called kinetic transformations [20]. Let us consider the non-biochemical system (2). The first step in bio-transforming (2) is to translate relevant trajectories5 x , ¯ y ) into the nonnegative quadrant. However, since R α ( X ) is a black-box network, i.e. f ( x ; α )is unknown and unalterable, only ¯ y can be translated; to this end, we define a new variable y ≡ (¯ y + T ), with translation T >
0, under which (2) becomesd x d t = f ( x ; α ) + γ y − γ T, x ∗ = β γ , d y d t = β − γ x , y ∗ = − γ − f (cid:18) β γ ; α (cid:19) + T. (5)Terms ( − γ x ) and ( − γ T ), called cross-negative terms [20], do not correspond to biochemicalreactions and, therefore, must be eliminated. Let us note that cross-negative terms also play acentral role in the questions of existence of other fundamental phenomena in biochemistry, suchas oscillations, multistability and chaos [20, 26]. Term ( − γ x ) can be eliminated with the so-called hyperbolic kinetic transformation, presented in Appendix D, which involves introducing anadditional controlling species Y and expanding system (5) intod x d t = f ( x ; α ) + γ y − γ T, x ∗ = β γ , d y d t = β − β y y , y ∗ = − γ − f (cid:18) β γ ; α (cid:19) + T, d y d t = γ x − β y y , y ∗ = β β ( y ∗ ) − . (6)Note that (5) and (6) have identical equilibria (time-independent solutions) for the species X and Y , and that the equilibria for the species Y and Y have a hyperbolic relationship; furthermore,provided β is sufficiently large, time-dependent solutions of (5) and (6) are close as well, seeAppendix D. On the other hand, cross-negative term ( − γ T ) can be eliminated via multiplicationwith x and any other desired factor, without influencing the x -equilibrium, which is determinedsolely by the RREs for y and y , and which we want to preserve. One option is to simply map( − γ T ) to ( − γ T x ), and take T large enough to ensure that the y ∗ -equilibrium is positive; however,this approach requires the knowledge of f ( x ; α ). A more suitable approach is to map ( − γ T ) to( − γ T x y ), under which, defining γ ≡ γ T , one obtainsd x d t = f ( x ; α ) + γ y − γ x y , x ∗ = β γ , d y d t = β − β y y , y ∗ ) + (cid:20) γ − f (cid:18) β γ ; α (cid:19)(cid:21) y ∗ − (cid:18) γ γ γ β β (cid:19) , d y d t = γ x − β y y , y ∗ = β β ( y ∗ ) − . (7)The quadratic equation for y ∗ from (7) always has one positive solution; therefore, there alwaysexists an equilibrium with positive y - and y -components.In what follows, we consider input networks R α ( X ) with at most two target species, { X , X } ,6nd we focus on controlling species X with the bimolecular controllers induced by (7), given by R β ( Y , Y ) : ∅ β −→ Y ,Y + Y β −→ ∅ , R γ ( Y ; X ) : X γ −→ X + Y , R + γ ( X i ; Y ) : Y γ −→ X i + Y , for some i ∈ { , } , R − γ ( X j ; Y ) : X j + Y γ −→ Y , for some j ∈ { , } . (8)In particular, the controller core R β ( Y , Y ) consists of a production of Y from a source, which wedenote by ∅ , and a bimolecular degradation of Y and Y . On the other hand, the controller interfaceconsists of the unimolecular reactions R γ ( Y ; X ), and R + γ ( X i ; Y ), that produce Y catalyticallyin X , and X i catalytically in Y , respectively, and the bimolecular reaction R − γ ( X j ; Y ) thatdegrades a target species X j catalytically in Y . We call reactions R + γ ( X i ; Y ) and R − γ ( X j , Y ) positive and negative interfacing, respectively. Furthermore, we say that positive (respectively,negative) interfacing is direct if i = 1 (respectively, if j = 1), i.e. if it is applied directly to thecontrolled species X ; otherwise, the interfacing is said to be indirect . In Figure 1, we displaycontroller (8) with direct positive and negative interfacing applied to an input network with asingle target species. Let us note that the AIFC from [22] is of the form (8), but it lacks negativeinterfacing R − γ ( X j ; Y ). In view of the derivation from this section, the AIFC is missing a keydesigning step, namely the translation from (5), indicating that a NEC may occur due to y ∗ - and y ∗ -equilibria being negative. Let us consider the unimolecular input network R α = R α ( X , X ), given by R α ( X , X ) : ∅ α − (cid:42)(cid:41) − α X , X α −→ X + X , X α −→ ∅ . (9)We interpret (9) as a simplified model for intracellular gene-expression. In this context, X isa degradable protein species that is produced via translation from a degradable mRNA species X , which is transcribed from a gene; species that are not explicitly modelled, such as genes,transcription factors and waste molecules, are denoted by ∅ . See also Figure 3(a) for a schematicrepresentation of network (9). The RREs of (9) have a unique globally stable equilibrium given by x ∗ = α α α α , x ∗ = α α . (10)The goal in this section is to control the equilibrium concentration of the protein species X at thedeterministic level, and its average copy-number at the stochastic level. To this end, we embeddifferent variants of the controller (8) into (9). Pure positive interfacing . Let us first consider controller (8) with only positive interfacing,i.e. the AIFC from [22]. We denote the controller by R + β,γ ≡ R β ∪ R γ ∪ R + γ and, for simplicity,7a) (b)Figure 3: Schematic representation of the gene-expression input network (9) . Panel (a) displays (9) with basal transcription rate α . Panel (b) display network (9) with tripled effective transcriptionrate, α , as a result of an activating transcription factor binding to the underlying gene promoter. assume that interfacing is direct: R β ( Y , Y ) : ∅ β −→ Y ,Y + Y β −→ ∅ , R γ ( Y ; X ) : X γ −→ X + Y , R + γ ( X ; Y ) : Y γ −→ X + Y . (11)The RREs for the output network (9) ∪ (11) are given byd x d t = α x − α x + γ y , d x d t = α − α x , d y d t = β − β y y , d y d t = γ x − β y y , (12)with the unique equilibrium x ∗ = β γ , x ∗ = α α , y ∗ = α γ (cid:18) β γ − α α α α (cid:19) , y ∗ = β β ( y ∗ ) − . (13)As anticipated in Section 3, the AIFC can lead to equilibria with negative y - and y -components.In particular, equation (13) implies that the nonnegative equilibrium is destroyed when β /γ <α α / ( α α ). Put it another way, using only positive interfacing, it is not possible to achievean output equilibrium below the input one. To determine the long term-behavior of (9) ∪ (11)when the nonnegative equilibrium ceases to exist, let us consider the linear combination of speciesconcentration ( α − x + α α − α − x + γ − ( y − y )) that, using (12), satisfiesdd t (cid:18) α x + α α α x + 1 γ ( y − y ) (cid:19) = − (cid:18) β γ − α α α α (cid:19) + γ α y ≥ − (cid:18) β γ − α α α α (cid:19) . (14)When β /γ < α α / ( α α ), equation (14) implies that a species concentration blows up for allnonnegative initial conditions, i.e. the output network displays a deterministic NEC; by applyingidentical argument to the first-moment equations, it follows that a stochastic NEC occurs as well.This result is summarized as a bifurcation diagram in Figure 4(a).To demonstrate these considerations, let us consider (9) with rate coefficients α = ( α , α , α , α )fixed so that the input equilibrium from (10) is given by x ∗ = 5. We choose the rate coefficientsof the AIFC (11) so that the target x -equilibrium from (13) is β /γ = 10; this setup is shown as8 (a) (e) (i) (b) (f) (j) (c) (g) (k) (d) (h) (l)Figure 4: Application of the IFC (8) on the input network (9) with rate coefficients ( α , α , α ) =(2 / , , α = α ( t ) which changes at t = 50 and leads to a catastrophicbifurcation . Panel (a) displays a bifurcation diagram for the output network (9) ∪ (11) , while panels (b)–(d) show the underlying deterministic and stochastic trajectories with α = 4 for t < and α = 12 for t ≥ , leading to a change in the parameter space indicated by the black arrow inpanel (a) . The control coefficients are fixed to ( β , β , γ , γ ) = (40 , , , . Analogous plots areshown in panels (e)–(h) for the output network (9) ∪ (15) with α = 12 for t < and α = 4 for t ≥ , and ( β , β , γ , γ ) = (40 , , , , and for the output network (9) ∪ (16) in panels (i)–(l) usingthe same coefficient values as in panels (a)–(d) , and with γ = 4 . ,
10) in Figure 4(a), which lies in the region where the output network displays anonnegative equilibrium. However, assume that at a future time, as a response to a signal fromthe environment, an activating transcription factor binds to the underlying gene promoter, triplingthe transcription rate of the input network, see Figure 3(b). Such a change would move the sys-tem from coordinate (5 ,
10) to (15 , ∪ (11), respectively, for the species X , Y and Y ; also shownas a dashed grey line in Figure 4(b) is the target equilibrium β /γ = 10. For time t <
50, whenthe input network operates as in Figure 3(a), and the output system is in the configuration (5 , t >
50, when the transcription hasincreased as in Figure 3(b), and the output system is at (15 ,
10) from Figure 4(a), control fails;even worse, the species Y blows up for all admissible initial conditions. Intuitively, this hazardousphenomenon occurs because, when the target equilibrium is below the input one, the best resultthat the AIFC can achieve, from the accuracy perspective, is to minimally increase X . Such atask is accomplished with y →
0, which, as a consequence of a hyperbolic relationship between y and y , enforces y → ∞ , which is a worst result from the stability perspective. Pure negative interfacing . Consider controller (8) with pure direct negative interfacing,denoted by R − β,γ ≡ R β ∪ R γ ∪ R − γ and given by R β ( Y , Y ) : ∅ β −→ Y ,Y + Y β −→ ∅ , R γ ( Y ; X ) : X γ −→ X + Y , R − γ ( X j ; Y ) : X + Y γ −→ Y . (15)By analogous arguments as with controller (11), one can prove that deterministic and stochasticNECs occur when β /γ > α α / ( α α ), i.e. controller (15) one can only decrease the input x -equilibrium from (10) to a lower output value. When control above the input equilibrium isattempted, species Y blows up. We display a bifurcation diagram and trajectories for the outputnetwork (9) ∪ (15) in Figure 4(e)–(h). Combined positive and negative interfacing . Let us now analyze controller (8) with bothpositive and negative interfacing directly applied to the protein species X , as suggested by thederivation in Section 3. This controller is denoted by R ± β,γ ≡ R β ∪ R γ ∪ R + γ ∪ R + γ and given by R β ( Y , Y ) : ∅ β −→ Y ,Y + Y β −→ ∅ , R γ ( Y ; X ) : X γ −→ X + Y , R + γ ( X ; Y ) : Y γ −→ X + Y , R − γ ( X j ; Y ) : X + Y γ −→ Y . (16)The RREs of the output network (9) ∪ (16) have two equilibria, given by x ∗ = β γ , x ∗ = α α , y ∗ = β β ( y ∗ ) − , (17)10here y ∗ satisfies ( y ∗ ) + α γ (cid:18) α α α α − β γ (cid:19) y ∗ − (cid:18) γ γ γ β β (cid:19) = 0 . (18)By design from Section 3, quadratic equation (18) always has exactly one positive equilibrium, sothat no NEC can occur with controller (16); we confirm this fact in Figure 4(i)–(l). Test network (9) demonstrates that controller (8) with only positive or negative interfacing does notgenerically ensure stability of the output network. In other words, the output network experiencesa NEC over a large range of parameters, as displayed in Figures 4(a) and (e). On the other hand,controller (8) with combined positive and negative interfacing, applied directly to the species ofinterest, induces no NEC, as shown in Figure 4(i). A special feature of unimolecular networks is thatdistinct species cannot influence each other negatively (via negative feedback loops). Consequently,to ensure existence of a nonnegative equilibrium, negative interfacing must be applied directly tothe target species whose dynamics is controlled, while positive interfacing can be applied directlyor indirectly. We now more formally state this result, which is proved in Appendix E. To aid thestatement of the theorem, consider the set of all unimolecular networks whose x -equilibrium iszero, x ∗ = 0, with the other target species X either present or deleted from the network. Thesenetworks form a negligibly small subset of general unimolecular networks, and are called degenerate;the set of all other unimolecular networks is said to be nondegenerate . Theorem 4.1.
Consider a nondegenerate unimolecular input network R α whose RRE s have anasymptotically stable equilibrium, the family of controllers R β,γ given by (8) , and the output network R α,β,γ = R α ∪R β,γ . Then, controller R ± β,γ with both positive and negative interfacing, with negativeinterfacing being direct, ensures that the output network R α,β,γ has a nonnegative equilibrium for allparameter values ( α , β , γ ) ∈ R a + b + c> . On the other hand, the other variants of the controller (8) donot generically ensure that the output network R α,β,γ has a nonnegative equilibrium; in the absenceof a nonnegative equilibrium, these controllers induce deterministic and stochastic blow-ups ( NEC s)for all nonnegative initial conditions.Proof.
See Appendix E.Let us note that the only variant of (8) that generically ensures a nonnegative equilibrium isalso the one which may be experimentally most challenging to implement. In particular, one mustgenerally implement the second-order reaction R − γ from (8) applied directly to the species whosedynamics is controlled.Let us also note that, in the degenerate case when the equilibrium of the target species fromthe input network is zero, x ∗ = 0, in addition to the controller (8) with both positive and negativeinterfacing, the AIFC also generically prevents NECs. However, as mentioned before, such inputnetworks are of limited relevance in applications. For example, when α = 0, the equilibriumof network (9) is zero and, consequently, the output equilibrium (13) is always nonnegative. Infact, identical gene-expression network without basal transcription, α = 0, has been used in [22]to demonstrate a desirable stochastic behavior of the AIFC. As we have shown in this section,when a more realistic gene-expression model is used, with α (cid:54) = 0, the AIFC can fail and induceboth deterministic and stochastic catastrophes. Similar degenerate input networks with a zero-equilibrium have also been used in [24, 25]. 11 Control of bimolecular input networks: Residual NEC
Most biochemical networks from applications are bimolecular, rather than unimolecular, limitingthe applicability of Theorem 4.1. In particular, Theorem 4.1 is rooted in a special property thatstable unimolecular networks controlled with (8) display: ensuring that the equilibrium of Y ispositive, y ∗ >
0, is sufficient to eliminate NECs, which is why positive and direct negative interfacinghave been sufficient to prevent NEC from occurring in Section 4. Let us now investigate if positivityof y ∗ is sufficient to eliminate NECs when bimolecular networks are controlled. To this end, considerthe input network R α ( X , X ), given by R α ( X , X ) : ∅ α −→ X , X α −→ X , X + X α −→ X , X α −→ ∅ , (19)where X is produced from a source and converted into a degradable species X via first- andsecond-order conversion reactions. We assume that X is a target species, while X is a residualspecies, i.e. X cannot be interfaced with a given controller. The RREs of (19) have a uniqueasymptotically stable equilibrium, given by x ∗ = α α α α + α α , x ∗ = α α . (20)(a) (b)(c) (d)Figure 5: Application of the IFC (16) on the input network (19) with rate coefficients ( α , α , α ) =(1 , , / α = 1 /
10 for t <
50, which changes to α = 3 /
10 for t ≥
50 and leads to acatastrophic bifurcation . Panels (a) – (d) display the deterministic and stochastic trajectories forthe output network (16) ∪ (19) , with control coefficients ( β , β , γ , γ , γ ) = (100 , , , , . Consider now controller R ± β,γ , given by (16), that always guarantees existence of a nonnegativeequilibrium for stable unimolecular input networks. Embedding (16) into (19), one can readily12how that the RREs of the resulting output network have two equilibria, one of which is nevernonnegative, while the other equilibrium satisfies x ∗ = β γ , x ∗ = α α β γ (cid:18) α α − β γ (cid:19) , y ∗ > , y ∗ > . (21)Therefore, output network (16) ∪ (19) has a nonnegative equilibrium provided β γ < α α . (22)Critically, even though y ∗ >
0, the equilibrium of the residual species X , that is not interfaced withthe controller, can become negative; furthermore, the combination of parameters α /α cannot besimply interpreted as a component of the input equilibrium (20).In Figure 5, we display the deterministic and stochastic trajectories for the output network (16) ∪ (19) over a time-interval such that condition (22) is satisfied for t <
50, and violated for t ≥
50 dueto a change in α . One can notice that the output network undergoes deterministic and stochasticNECs. Critically, not only does the controlling species Y blow up, but also the residual species X . Put it another way, controller (16) destabilizes the originally asymptotically stable inputnetwork (19). We call this hazardous phenomenon, arising from the fact that a residual specieshas no nonnegative equilibria, a residual NEC. Intuitively, when the concentration of the targetspecies X is increased beyond the upper bound from (22), residual species X , which influences X negatively, counteracts the positive action of the controlling species Y , resulting in a joint blow-up.We also display this phenomenon in context of intracellular control in Figure 2.Test network (20) demonstrates that, as opposed to unimolecular networks, positivity of y ∗ is not always sufficient to eliminate NECs when bimolecular networks are controlled, because theunderlying nonnegative equilibria of the residual species can cease to exist. In context of designingbiochemical networks [20], such a problem is overcome by translating not only the controllingspecies concentration y , but also the desired trajectories of the input (including residual) species x . However, as outlined in Section 3, in context of controlling molecular networks, the inputnetwork is unknown and unalterable, so that translations of x cannot be achieved and, therefore,residual NECs cannot be eliminated. When embedded into black-box biochemical networks, synthetic molecular controllers must achievea desired accuracy, while preserving stability. In this paper, we have demonstrated that the abilityof molecular IFCs to robustly achieve perfect accuracy comes at a severe stability cost. In particular,under these controllers, all nonnegative equilibria of a given input network can be lost and some ofthe underlying species abundances can blow up at both deterministic and stochastic levels. We callthis biochemically hazardous phenomenon a negative-equilibrium catastrophe (NEC). In context ofelectro-mehanical systems, analogous phenomenon, involving some variables reaching the boundaryof physically allowed values, is known as integrator windup [19].We have shown in Section 2 that, in contrast to linear electro-mechanical systems, unimolecularbiochemical systems cannot achieve integral control because of NECs. We have then constructed afamily of bimolecular IFCs (8) in Section 3 which, due to a careful designing procedure, affect inputnetworks both positively and negatively. In Section 4, we have shown that these dual-acting IFCs,whose negative interfacing is applied directly to the controlled biochemical species, are safeguardedagainst NECs when applied to any stable unimolecular network. However, in Section 5, we have13emonstrated that this safety profile is lost when bimolecular input networks are considered. Inparticular, when stable bimolecular networks are controlled, nonnegative equilibria for the residual species, that are not interfaced with the controller, can vanish, and the species blow up. This moresevere, destabilizing, form of NEC is called residual
NEC, and it cannot be eliminated without adetailed knowledge about the structure and kinetics of the input network.Most realistic biochemical networks from applications are bimolecular, and contain a largenumber of coupled biochemical species, the vast majority of which are residual (not interfaceablewith a given controller), see also Figure 1. When such networks are controlled with molecular IFCs,it takes only one of the many residual species to not display a nonnegative equilibrium for the controlto fail and a catastrophic event to unfold. This phenomenon has already been demonstrated withnetwork (19), which contains only one bimolecular reaction and only one residual species. Therefore,residual NECs place a fundamental limit to applicability of molecular IFCs in synthetic biology.In particular, without a detailed knowledge of the controlled systems, no molecular IFC can besafely applied to bimolecular, and hence the vast majority of intracellular, networks, without therisk of causing a catastrophic blow-up of molecular abundances. Furthermore, even if equippedwith a detailed knowledge about an input network, determining the parameter regimes that ensureexistence of a nonnegative equilibrium can be a formidable task involving solving a large numberof nonlinear equations, and experimentally fine-tuning control coefficients with a high precision.
Acknowledgements
This work was supported by the EPSRC grant EP/P02596X/1. Guy-Bart Stan acknowledges sup-port from his UK Royal Academy of Engineering Chair in Emerging Technologies for EngineeringBiology (RAE CiET1819 \ A Appendix: Background
Notation . Given sets A and A , their union, intersection, and difference are denoted by A ∪ A , A ∩ A , and A \ A , respectively. The largest element of a set of numbers A is denoted bymax A . The empty set is denoted by ∅ . Set Z is the space of integer numbers, Z ≥ the space ofnonnegative integer numbers, and Z > the space of positive integer numbers. Similarly, R is thespace of real numbers, R ≥ the space of nonnegative real numbers, and R > the space of positivereal numbers. Euclidean column vectors are denoted in boldface, x = ( x , x , . . . , x N ) (cid:62) ∈ R N = R N × , where · (cid:62) denotes the transpose operator. The i -standard basis vector is denoted by e i ≡ ( δ i, , δ i, , . . . , δ i,N ) (cid:62) ∈ R N , where δ i,j = 0 if i (cid:54) = j , and δ i,j = 1 if i = j . The zero-vector is denoted by ≡ (0 , , . . . , (cid:62) ∈ R N , and we let ( ≡ (cid:80) Ni =1 e i ) ∈ R N . Given two Euclidean vectors x , y ∈ R N ,their inner-product is denoted by (cid:104) x , y (cid:105) ≡ (cid:80) Ni =1 x i y i . Abusing the notation slightly, given twosequences u, v : Z ≥ → R , their inner-product is also denoted by (cid:104) u ( x ) , v ( x ) (cid:105) ≡ (cid:80) x ∈ Z ≥ u ( x ) v ( x );we let (cid:107) u ( x ) (cid:107) l ≡ (cid:80) x ∈ Z ≥ | u ( x ) | denote the l -norm of u ( x ). Given a matrix A ∈ R N × N , with( i, j )-element α i,j ∈ R , we denote the i -row and j -column of A by α i, · ≡ ( α i, , α i, , . . . , α i,N ) ∈ R × N ≥ , and α · ,j ≡ ( α ,j , α ,j , . . . , α N,j ) (cid:62) ∈ R N × ≥ , respectively. The identity matrix is denoted by I ≡ ( e , e , . . . , e N ) ∈ R N × N , the zero matrix by 0 ≡ ( , , . . . , ) ∈ R N × N , and diagonal matricesare denoted by diag( α , , α , , . . . , α N,N ) ≡ ( α , e , α , e , . . . , α N,N e N ) ∈ R N × N . Determinant ofa matrix A is denoted by | A | . 14 .1 Biochemical reaction networks Let R α = R α ( X ) be a reaction network describing interactions, under mass-action kinetics, between N biochemical species X = { X , X , . . . , X N } in a well-mixed unit-volume reactor [10], as specifiedby the following a reactions: R α ( X ) : N (cid:88) l =1 ν j,l X l α j −→ N (cid:88) l =1 ¯ ν j,l X l , j ∈ A = { , , . . . , a } . (23)Here, α = ( α , α , . . . , α A ) ∈ R a> are the positive rate coefficients of the reactions from R α . Non-negative vectors ν j, · = ( ν j, , ν j, , . . . , ν j,N ) (cid:62) ∈ Z N ≥ and ¯ ν j, · = (¯ ν j, , ¯ ν j, , . . . , ¯ ν j,N ) (cid:62) ∈ Z N ≥ are the reactant and product stoichiometric coefficients of the j -reaction, respectively; if ν j, · = (respec-tivelty, ¯ ν j, · = ), then the reactant (respectively, product) of the j -reaction is the null-species , de-noted by ∅ , representing species that are not explicitly modelled. When convenient, we denote twoirreversible reactions ( (cid:80) Nl =1 ν i,l X l α i −→ (cid:80) Nl =1 ¯ ν i,l X l ) ∈ R α and ( (cid:80) Nl =1 ¯ ν i,l X l α j −→ (cid:80) Nl =1 ν i,l X l ) ∈ R α jointly as the single reversible reaction ( (cid:80) Nl =1 ν i,l X l α i − (cid:42)(cid:41) − α j (cid:80) Nl =1 ¯ ν i,l X l ) ∈ R α . Species X i is a catalyst in the j -reaction from R α if ν j,i = ¯ ν j,i (cid:54) = 0; if X i is a catalyst in all of the reaction from R α , then wewrite R α = R α ( X \ X i ; X i ). The order of the j-reaction from network R α is given by (cid:104) , ν j, · (cid:105) ∈ Z ≥ .The order of reaction network R α is given by max {(cid:104) , ν j, · (cid:105)| j ∈ A} ; first-order (respectively, second-order) reaction networks are also said to be unimolecular (respectively, bimolecular ).Given a class of biochemical reaction networks, parametrized by the underlying rate coefficients,it may be of interest if a given property is likely to be true when all admissible values of the ratecoefficients are considered, which motivates the following definition. In what follows, we implicitlyuse Lebesgue measure for sets. Definition A.1 ( Genericity ) . Consider a mass-action reaction network R α parametrized by therate coefficients α ∈ S α , where S α ⊂ R a> is a nonempty open set. Assume S α is partitioned accordingto S α = Ω α ∪ ω α , with Ω α ∩ ω α = ∅ , where ω α is a set of measure zero. A property is said to be generic for the set S α and network R α if it holds for all α ∈ Ω α and fails to hold for all α ∈ ω α .Example . Empy set ω α = ∅ , and set ω α = { α , α , . . . } , containing finitely or countably infinitelymany points, have measure zero. All the points α where a non-trivial polynomial P ( α ) vanishes isalso a set of measure zero [28]; for example, given a nonzero matrix A = A ( α ), with ( i, j )-element α i,j , the set of all points α such that | A ( α ) | = 0 has zero measure. A.2 Dynamical models of reaction networks
In what follows, we present deterministic and stochastic models of mass-action reaction networks,and provide definitions in context of blow-ups.
A.2.1 Deterministic model
Let x ( t ; α ) = ( x ( t ; α ) , x ( t ; α ) , . . . , x N ( t ; α )) (cid:62) ∈ R N ≥ be a concentration vector at time t ∈ R ≥ for the species X = { X , X , . . . , X N } from the network R α , given by (23). A deterministic modelof the reaction network R α describes the time-evolution of x = x ( t ; α ) as a system of first-orderordinary differential equations (ODEs), called the reaction-rate equations (RREs) [10, 9], given byd x d t = K ( x ; α ) = (cid:88) j ∈A α j ∆ x j, · x ν j, · , (24)15here ∆ x j, · ≡ (¯ ν j, · − ν j, · ) ∈ Z N is the reaction vector of the j -reaction, and x ν j, · ≡ (cid:81) Ni =1 x ν j,i i with 0 ≡
1. Function K ( · ; α ) : R N → R N , called a kinetic function (see also Appendix D), is apolynomial of degree m = max {(cid:104) , ν j, · (cid:105)| j ∈ A} in x , which we denote by K ( x ; α ) ∈ P m ( R N ; R N ).Vector x ∗ = x ∗ ( α ) < ∞ is called an equilibrium of the RREs (24) if K ( x ∗ ; α ) = .We now present an important property of the RREs [20]. Theorem A.1.
The nonnegative orthant R N ≥ is an invariant set for the ODE s (24) .Proof. See [20].In this paper, unimolecular networks are of interest and, to this end, we introduce the followingdefinition.
Definition A.2 ( Cross-nonnegative matrix ) . A matrix A ∈ R N × N with nonnegative off-diagonal elements, α i,j ≥ for all i, j ∈ { , , . . . , N } such that i (cid:54) = j , is said to be cross-nonnegative .Remark . Cross-nonnegative matrices are known as negative Z , quasi-positive, essentially nonneg-ative, and Metzler matrices in the literature [29, 30].In context of unimolecular reaction networks, Theorem A.1 implies the following corollary. Corollary A.1.
For unimolecular reaction networks, the kinetic function from the
RRE s (24) isgiven by K ( x ; α ) = ( α · , + A x ) , where α · , = ( α , , α , , . . . , α N, ) ∈ R N ≥ is a nonnegative vector,while A ∈ R N × N is a cross-nonnegative matrix. A linear ODE system d x / d t = ( α · , + A x ) is said to be asymptotically stable (also simplyreferred to as stable in this paper) if all the eigenvalues of A have negative real parts. A.2.2 Stochastic model
Let X ( t ; α ) = ( X ( t ; α ) , X ( t ; α ) , . . . , X N ( t ; α )) (cid:62) ∈ Z N ≥ be a copy-number vector at time t ∈ R ≥ for the species X = { X , X , . . . , X N } from the network R α , given by (23); abusing the nota-tion slightly, we denote the points in the state-space for X ( t ; α ) using the same symbol as theconcentration vector from (24), i.e. by x = ( x , x , . . . , x N ) (cid:62) ∈ Z N ≥ . A stochastic model of the re-action network R α describes the time-evolution of X = X ( t ; α ) as a continuous-time discrete-spaceMarkov chain [31] characterized via a difference-operator L α , called the generator , given by [32] L α u ( x ) = A (cid:88) j =1 α j x ν j, · ( E +∆ x j, · x − u ( x ) , (25)where u : Z N ≥ → R belongs to a suitable function space. Here, ∆ x j, · = (¯ ν j, · − ν j, · ) ∈ Z N is the j -reaction vector, while x ν j, · = (cid:81) Ni =1 x ν j,i i , with x ν j,i i = x i ( x i − . . . ( x i − ν j,i −
1) and x ≡ x ∈ Z ≥ . Furthermore, E +∆ x j, · x = (cid:81) Ni =1 E +∆ x j,i x i is a step-operator such that E +∆ x j, · x u ( x ) = u ( x + ∆ x j, · ).Let p ( · , t ; α ) : Z N ≥ → [0 ,
1] be the probability-mass function (PMF) at time t of the Markovchain with generator (25), and let f x : Z N ≥ → R be a suitable function of the species copy-numbers.We let E f x = E f x ( t ; α ) ≡ (cid:104) f x ( x ) , p ( x , t ; α ) (cid:105) denote the expectation of f x at time t with respect to p ( x , t ; α ). In this paper, we focus on the average species copy-numbers, i.e. on the first-moment16ector E X = E X ( t ; α ) = ( E X ( t ; α ) , E X ( t ; α ) , . . . , E X N ( t ; α )) (cid:62) ∈ R N ≥ , which evolves in timeaccording to the ODEs [9, 32] d E X d t = E [ L α X ] = (cid:88) j ∈A α j ∆ x j, · E X ν j, · . (26)In the special case of unimolecular reaction networks, the first-moment equations (26) and theRREs (24) are formally equivalent; more generally, the less-detailed deterministic and the more-detailed stochastic models match in the thermodynamic limit [33]. A.2.3 Blow-up and negative-equilibrium catastrophe
In this paper, we focus on the circumstances when some of the species abundance experiences anunbounded growth, and introduce the following definition for this purpose.
Definition A.3 ( Blow-up ) . Reaction network R α ( X ) is said to blow up deterministically fora given initial condition if lim t →∞ x i ( t ; α ) = ∞ for some i ∈ { , , . . . , N } , where the speciesconcentration x ( t ; α ) ∈ R N ≥ satisfies (24) ; R α ( X ) is said to blow up stochastically for a giveninitial condition if lim t →∞ E X i ( t ; α ) = ∞ for some i ∈ { , , . . . , N } , where the first-moment ofthe species copy-numbers E X ( t ; α ) ∈ R N ≥ satisfies (26) . Nonnegative ODEs, such as equations (24) and (26), need not have a nonnegative time-independentsolution, i.e. the dynamics can be confined to an unbounded invariant set devoid of any equilibria.In this context, we introduce the following definition.
Definition A.4 ( Negative-equilibrium catastrophes (NECs) ) . Reaction network R α ( X ) issaid to display a deterministic (respectively, a stochastic) negative-equilibrium catastrophe ( NEC )if, for all α ∈ R a> such that the RRE s have no nonnegative equilibria, R α ( X ) blows up determin-istically (respectively, stochastically) for some nonnegative initial conditions.Remark . A nonnegative equilibrium can cease to exist by attaining a negative component, becomingcomplex, or vanishing all together. B Appendix: Stochastic biochemical control
In this section, we formulate the problem of achieving biochemical control over a given reactionnetwork, starting with the following definition.
Definition B.1 ( Black, grey and white box ) . Network R α ( X ) (cid:54) = ∅ with unknown (respectively,only partially known) structure and dynamics is called a black-box (respectively, grey-box ) network; R α ( X ) (cid:54) = ∅ is called a white-box network if its structure and dynamics are completely known. Given a black- or grey-box input (uncontrolled) reaction network R α = R α ( X ), the objec-tive of biochemical control is to design a controller network R β,γ in order to ensure that thedynamics of desired input species X is suitably controlled in the resulting output (controlled) net-work R α,β,γ ≡ R α ∪ R β,γ . To this end, we partition the input species into X = X τ ∪ X ρ , where X τ = { X , X , . . . , X n } , with 1 ≤ n ≤ N , are the target species that can be interfaced with agiven controller, while X ρ = X \ X τ = { X n +1 , X n +2 , . . . , X N } are the residual species that can-not be interfaced with the controller. The controller can be decomposed into two sub-networks, R β,γ = R β,γ ( X τ , Y ) = R β ( Y ) ∪ R γ ( X τ , Y ), where R β = R β ( Y ), called the core , contains all the17eactions that involve only the controlling species Y = { Y , Y , . . . , Y M } , while R γ = R γ ( X τ , Y ),called the interface , contains all of the remaining reactions, involving both X τ and Y . Put moresimply, the core R β describes internal dynamics of the controlling species, while the interface R γ describes how the target and controlling species interact.In this paper, we focus on controlling average copy-number of a single target species; see [34] fora more general multi-species control of the full PMF of both target and residual species. To this end,we denote the rate coefficients from the sub-networks R α , R β and R γ by α ∈ R a> , β ∈ R b> and γ ∈ R c> , respectively. We also let E X ( t ; · , · ; α ) : R b> × R c> → R be the first-moment of the target species X ∈ X τ , which is a function of the control parameters ( β , γ ) for every time t , input coefficient α and every initial condition. More precisely, E X = E X ( t ; β , γ ; α ) = (cid:104) x , p ( x , y , t ; β , γ ; α ) (cid:105) , where p ( x , y , t ; β , γ ; α ) is the time-dependent PMF of the output network. In what follows, we denotethe gradient operator with respect to x = ( x , x , . . . , x N ) by ∇ x ≡ ( ∂/∂ x , ∂/∂ x , . . . , ∂/∂ x N ). Definition B.2 ( Control ) . Consider a black-box input network R α ( X ) , and the correspondingoutput network R α,β,γ ( X , Y ) = R α ( X ) ∪ R β,γ ( X τ , Y ) . Assume we are given a nonempty open set S α ⊂ R a> , stability index p ∈ Z > , and a target value ¯ x ∈ R > together with a tolerance ε ∈ R ≥ .Then, the first-moment E X = E X ( t ; β , γ ; α ) is said to be controlled in the long-run if thereexists a nonempty open set S β,γ ⊂ R b + c> such that the following two conditions are satisfied for allinitial conditions ( X (0) , Y (0)) ∈ Z N + M ≥ and rate coefficients ( α , β , γ ) ∈ S α × S β,γ : (C.I) Stability : { lim t →∞ E X pi ( t ; β , γ ; α ) < ∞} Ni =1 and { lim t →∞ E Y pi ( t ; β , γ ; α ) < ∞} Mi =1 . (C.II) Accuracy : lim t →∞ (cid:12)(cid:12)(cid:12) E X ( t ; β , γ ; α ) − ¯ x (cid:12)(cid:12)(cid:12) ≤ ε , where lim t →∞ ∇ β , γ E X ( β , γ ; α ) (cid:54) = . Condition (C.I) requires boundedness of the long-time moments, up to order p >
1, of all ofthe species from the output network R α,β,γ ( X , Y ). Minimally, one requires that the long-timefirst-moments are bounded, i.e. that the controller does not trigger a stochastic blow-up (seeDefinition A.3). The larger p ≥ X , which is required to depend on at least one control parameter, is sufficientlyclose to the target value. Let us remark that (C.I)–(C.II) must hold within neighborhoods of theunderlying rate coefficient values, reflecting the fact that measurement and fine-tuning of ratecoefficients is not error-free. For the same reason, we have put forward a more relaxed accuracycriterion by allowing nonzero tolerance in condition (C.II). B.1 Robust control
Condition (C.II) from Definition B.2 may be challenging to achieve due to the fact that the long-time first-moment E X generally depends on the initial conditions, and on the input coefficients α , which motivates the following definition. Definition B.3 ( Robustness ) . Consider a black-box input network R α ( X ) , and the correspond-ing output network R α,β,γ ( X , Y ) = R α ( X ) ∪ R β,γ ( X τ , Y ) . Assume conditions (C.I)–(C.II) from Definition B.2 are satisfied. Then, network R β,γ ( X τ , Y ) is a robust controller of the first-moment E X in the long-run if the following two conditions are also satisfied: (R.I) Robustness to initial conditions . There exists a unique stationary
PMF p ( x , y ; β , γ ; α ) such that lim t →∞ (cid:107) p ( x , y , t ; β , γ ; α ) − p ( x , y ; β , γ ; α ) (cid:107) l = 0 for all ( X (0) , Y (0)) ∈ Z N + M ≥ and ( α , β , γ ) ∈ S α × S β,γ . Robustness to input coefficients . The stationary first-moment, given by E X ∗ ( β , γ ; α ) ≡ (cid:104) x , p ( x , y ; β , γ ; α ) > , satisfies ∇ α E X ∗ ( β , γ ; α ) = for all ( α , β , γ ) ∈ S α × S β,γ .Remark . Robust controllers are also called integral-feedback controllers [19]. Remark . Conditions (C.II) and (R.II) demand that the first-moment of X is nondegenerate withrespect to the control parameters, ∇ β , γ E X ∗ ( β , γ ; α ) (cid:54) = , and that it is degenerate with respectto the input parameters, ∇ α E X ∗ ( β , γ ; α ) = , respectively. As we prove in Lemma C.1 below,the degeneracy condition enforces singularity of an appropriate kinetic matrix. C Appendix: Nonexistence of unimolecular integral-feedback con-trollers
Let R α ( X ) be a black-box input network with a desired target species X ∈ X τ , and let R β,γ ( X τ , Y ) = R β ( Y ) ∪ R γ ( X τ , Y ) be a unimolecular network. The first-moment equations for the output network R α,β,γ ( X , Y ) = R α ( X ) ∪ R β,γ ( X τ , Y ) can be written in the following form:d E X ρ d t = E [ L α X ρ ] , d E X τ d t = E [ L α X τ ] + C , E X τ + C , E Y , d E Y d t = β · , + C , E X τ + ¯ C , E Y , (27)where L α is the (unknown) generator of the input network, β · , = ( β , , β , , . . . , β M, ) ∈ R M ≥ isinduced by the core network R β ( Y ), while cross-nonnegative matrix C , ∈ R n × n and nonnegativematrices C , ∈ R n × M ≥ and C , ∈ R M × n ≥ are induced by the interfacing network R γ ( X τ , Y ). Matrix¯ C , ∈ R M × M ≥ can be written as ¯ C , = ( B + C , ), with cross-nonnegative matrices B and C , being induced by R β ( Y ) and R γ ( X τ , Y ), respectively. Note that ¯ C , encodes all of the reactionsthat change the copy-numbers of the controlling species Y , either in X τ -independent (via matrix B ) or X τ -dependent manner (via matrix C , ). Let us also note that, by definition, β ∈ R b> contains nonzero elements from vector β · , , and absolute value of nonzero elements from matrix B ;similarly, γ ∈ R c> contains absolute values of nonzero elements from C , , C , , C , and C , . Wenow present conditions on the matrices C , = ( γ , · , , γ , · , , . . . , γ , · ,n ) and ¯ C , that are necessary forunimolecular networks to exert robust (integral-feedback) control. Lemma C.1.
Assume a unimolecular reaction network R β,γ ( X τ , Y ) is a robust controller for ablack-box network R α ( X ) . Then, ¯ C , from (27) is a singular matrix all eigenvalues of which havenonpositive real parts for all ( β , γ ) ∈ S β,γ . Furthermore, the first column of matrix C , is nonzero.Proof. Assume ¯ C , has an eigenvalue with a positive real part for some ( β , γ ) ∈ S β,γ ; it thenfollows from (27) that lim t →∞ E Y i = ∞ for some i ∈ { , , . . . , M } . Therefore, condition (C.I) fromDefinition B.3 does not hold for all ( β , γ ) ∈ S β,γ .Assume ¯ C , is nonsingular for some ( β , γ ) ∈ S β,γ . Let ( E X ∗ τ , E X ∗ ρ , E Y ∗ ) be the unique sta-tionary first-moments, obtained by setting to zero the left-hand side in (27), and in the ODEs forhigher-order moments. Then, one can eliminate E Y ∗ via E Y ∗ = − ( ¯ C , ) − ( β · , + C , E X ∗ τ ). Theremaining system of equations for ( E X ∗ τ , E X ∗ ρ ), and hence E X ∗ , depends on α via the unknowngenerator L α . Therefore, condition (R.II) from Definition B.3 does not hold for all ( β , γ ) ∈ S β,γ .Assume the first column of C , is zero, γ , · , = . By the Fredholm alternative theorem, anecessary condition for the stationary first-moment E Y ∗ to exist is then given by ( (cid:104) w , β · , (cid:105) +19 ni =2 (cid:104) w , γ , · ,i (cid:105) E X ∗ i ) = 0 for every w ∈ R M such that ( ¯ C , ) (cid:62) w = . The stationary average E X ∗ isnot uniquely determined by these constraints and, hence, depends on α . Therefore, condition (R.II)from Definition B.3 does not hold. Theorem C.1.
There does not exist a unimolecular integral-feedback controller.Proof.
Suppose, for contradiction, that R β,γ ( X τ , Y ) is a robust (integral-feedback) controller, i.e.that the corresponding output network R α,β,γ ( X , Y ) = R α ( X ) ∪R β,γ ( X τ , Y ) satisfies Definition B.3.Lemma (C.1) then implies that ( ¯ C , ) (cid:62) is a singular cross-nonnegative matrix all eigenvalues ofwhich have nonpositive real parts; in what follows, we assume that ( ¯ C , ) (cid:62) is irreducible. It thenfollows from [29, Theorem 5.6] that there exists a positive vector w ∈ R N> such that ( ¯ C , ) (cid:62) w = ,so that d / d t (cid:104) w , E Y (cid:105) = ( (cid:104) w , β · , (cid:105) + (cid:80) ni =1 (cid:104) w , γ , · ,i (cid:105) E X i ). Using the fact that, by Lemma (C.1), γ , · , (cid:54) = , and that, by Definition B.3, lim t →∞ E X ( t ; β , γ ; α ) = E X ∗ ( β , γ ) >
0, it follows thatd / d t (cid:104) w , E Y (cid:105) > t →∞ E Y i = ∞ for some i ∈ { , , . . . , M } , implyingthat condition (C.I) from Definition B.2 does not hold. Hence, R β,γ ( X τ , Y ) is not a robust controller.If ¯ C , is reducible, then analogous argument can be applied to each of the underlying irreduciblecomponents, implying the statement of the theorem. Due to linearity of the controller, the sameargument implies that a deterministic integral-feedback controller does not exist.The proof of Theorem C.1 shows that a unimolecular integral-feedback controller does not existbecause the underlying reaction network R β,γ ( X τ , Y ) triggers a deterministic and stochastic NEC(see also Definition A.4 in Section A.2.3) for all nonnegative initial conditions, violating the stabilitycondition (C.I) from Definition B.2. For example, taking C , = ( γ , · , , , . . . , ), it follows that thedeterministic equilibrium of the target species is x ∗ = −(cid:104) w , β · , (cid:105) / (cid:104) w , γ , · , (cid:105) < D Appendix: Hyperbolic kinetic transformation
In this section, we briefly discuss how to map arbitrary polynomial ODEs into dynamically similarmass-action RREs [20]. To this end, let us note that every component of every polynomial function P ( x ; α ) = ( P ( x ; α ) , P ( x ; α ) , . . . , P N ( x ; α )) (cid:62) ∈ P m ( R N ; R N ), where α ∈ R a> , can be written as P i ( x ; α ) = (cid:88) j ∈A + i α j x ν j, · − (cid:88) j ∈A − i, K α j x ν j, · − (cid:88) j ∈A − i, N α j x ν j, · , ∀ i ∈ { , , . . . , N } . (28)Here, A + i and A − i ≡ ( A − i, K ∪ A − i, N ) are the indices of all of the distinct positive and negative termsin P i ( x ; α ), respectively, where A − i, K = { j ∈ A − i | ν i,j (cid:54) = 0 } and A − i, N = { j ∈ A − i | ν i,j = 0 } . Definition D.1 ( Cross-negative terms; kinetic functions ) . Consider a polynomial func-tion P ( x ; α ) ∈ P m ( R N ; R N ) with α ∈ R a> , whose i -component is given by (28) . Monomials {− α j x ν j, · } j ∈A − i, N are called cross-negative terms. Polynomial functions without cross-negativeterms, A − i, N = ∅ , are called kinetic functions and denoted by K ( x ; α ) ; the set of all kinetic func-tions of degree at most m is denoted by P K m ( R N ; R N ) . Polynomial functions with a cross-negativeterm, A − i, N (cid:54) = ∅ , are called non-kinetic functions and denoted by N ( x ; α ) ; the set of all non-kineticfunctions of degree at most m is denoted by P N m ( R N ; R N ) .Remark . Definition D.1 captures the fact that biochemical reactions cannot consume a species whenthe concentration of that species is zero, e.g. see network (1); Theorem A.1 is a direct consequenceof the absence of cross-negative terms in kinetic functions [20].20n order to map an arbitrary input non-kinetic function N ( x ; α ) ∈ P N m ( R N ; R N ) into an outputkinetic one K (¯ x ; ¯ α ) ∈ P K ¯ m ( R ¯ N ; R ¯ N ), while preserving desired dynamical features over suitable time-intervals, a number of so-called kinetic transformations have been developed in [20]. Such mappingsinvolve a dimension expansion (i.e. an introduction of additional biochemical species, ¯ N ≥ N ) oran increase in non-linearity ( ¯ m ≥ m ). We now present one such kinetic transformation. Definition D.2 ( Hyberbolic transformation ) . Consider a system of polynomial
ODE s givenby d x i d t = (cid:88) j ∈A + i α j x ν j, · − (cid:88) j ∈A − i, K α j x ν j, · , x i (0) = x i ≥ , for i ∈ { , , . . . , n } , d x i d t = (cid:88) j ∈A + i α j x ν j, · − (cid:88) j ∈A − i, K α j x ν j, · − (cid:88) j ∈A − i, N α j x ν j, · , x i (0) = x i > , for i ∈ { n + 1 , n + 2 , . . . , N } , (29) with the index sets as defined in (28) . Consider also the RRE s given by d x i d t = (cid:88) j ∈A + i α j x ν j, · − (cid:88) j ∈A − i, K α j x ν j, · , x i (0) = x i , for i ∈ { , , . . . , n } , d x i d t = (cid:88) j ∈A + i α j x ν j, · − (cid:88) j ∈A − i, K α j x ν j, · − ε x i ¯ x i , x i (0) = x i , for i ∈ { n + 1 , n + 2 , . . . , N } , d¯ x i d t = (cid:88) j ∈A − i, N α j x ν j, · − ε x i ¯ x i , ¯ x i (0) > , for i ∈ { n + 1 , n + 2 , . . . , N } , (30) with a parameter ε ∈ R > . Kinetic transformation Ψ εH : P m ( R N ; R N ) → P K ¯ m ( R (2 N − n ) ; R (2 N − n ) ) ,mapping the right-hand side of ODE s (29) to the right-hand side of (30) , with ¯ m ≤ m + 2 , is calleda hyberbolic transformation . By comparing equilibria of (29) and (30), the following proposition is established.
Proposition D.1.
Let x ∗ = ( x ∗ , x ∗ , . . . , x ∗ N ) ∈ R N be an equilibrium of the ODE system (29) with x ∗ i (cid:54) = 0 for all i ∈ { n + 1 , n + 2 , . . . , N } . Then, ( x ∗ , ¯x ∗ ) ∈ R N − n is an equilibrium of the ODE system (30) for some ¯x ∗ ∈ R N − n . Let us note that, while the x -equilibria are preserved by the hyperbolic kinetic transformation,their properties, such as stability, are not necessarily preserved; a stronger dynamical preservationis ensured by taking ε sufficiently small. Theorem D.1.
Solutions of (29) , with ( x n +1 , x n +2 , . . . , x N ) ∈ R N − n> , are asymptotically equivalentto the solutions of (30) in the limit ε → .Proof. Let us introduce a change of coordinates y i = ( x i − ¯ x i ) for i ∈ { n +1 , n +2 , . . . , N } , leading toa regular singularly perturbed system in the variables ( x , x , . . . , x n , y n +1 , y n +2 , . . . , y N , ¯ x n +1 , ¯ x n +2 ,. . . , ¯ x N ). The adjoined sub-system has an isolated equilibrium (¯ x ∗ n +1 , ¯ x ∗ n +2 , . . . , ¯ x ∗ N ) = which,when substituted into the degenerate sub-system in the variables ( x , x , . . . , x n ,y n +1 , y n +2 , . . . , y N ),leads to (29). If ( x n +1 , x n +2 , . . . , x N ) ∈ R N − n> , then the trivial adjoined equilibrium is stable, andTikhonov’s theorem [35] implies the statement of Theorem D.1.21 Appendix: Robust control of unimolecular input networks
In this section, we analyze performance of the class of second-order integral-feedback controllers R β,γ ( { X , X } , { Y , Y } ) = R β ( Y , Y ) ∪ R γ ( Y ; X ) ∪ R + γ ( X i ; Y ) ∪ R − γ ( X j ; Y ), given by (8), byembedding them into the class of unimolecular input networks whose RREs, given byd x d t = α · , + A x , (31)have an asymptotically stable equilibrium. We also assume the input network has two target species X τ = { X , X } , and the goal is to control the concentration/first-moment of the target species X .In what follows, we let A ( i ,i ,...,i n ) , ( j ,j ,...,j m ) ∈ R ( N − n ) × ( N − m ) denote the sub-matrix obtainedby removing from A ∈ R N × N the rows { i , i , . . . , i n } ⊂ { , . . . , N } and columns { j , j , . . . , j m }⊂ { , . . . , N } . Furthermore, we let A α · ,j → x ∈ R N × N denote the matrix obtained by replacing the j -column of A by a vector x ∈ R N . Theorem E.1. (Deterministic equilibria)
Let R α = R α ( X ) , with species X = { X , X , . . . , X N } ,be the class of unimolecular input networks whose RRE s (31) have an asymptotically stable equilib-rium. Let R β,γ = R β,γ ( { X , X } , { Y , Y } ) be the controller given by (8) . Then, the output network R α,β,γ = R α ∪ R β,γ satisfies the following properties: (i) Pure positive interfacing . Let R − γ = ∅ . If positive interfacing is direct, R + γ = R + γ ( X , Y ) (cid:54) = ∅ , then there exists a nonnegative equilibrium if and only if β γ > | A α · , →− α · , || A | ; the same istrue if positive interfacing is indirect, R + γ = R + γ ( X , Y ) (cid:54) = ∅ , and if | A , | (cid:54) = 0 . (ii) Pure negative interfacing . Let R + γ = ∅ . If negative interfacing is direct, R − γ = R − γ ( X , Y ) (cid:54) = ∅ , then there exists a nonnegative equilibrium if and only if β γ < | A →− α · , || A | . If negative inter-facing is indirect, R − γ = R − γ ( X , Y ) (cid:54) = ∅ , and if | A , | (cid:54) = 0 , then there exists a nonnegativeequilibrium if and only if | ( A α · , →− α · , ) , || A , | < β γ < | A α · , →− α · , || A | . (iii) Combined indirect negative interfacing . Let R − γ = R − γ ( X , Y ) (cid:54) = ∅ . If positive interfac-ing is direct, R + γ = R + γ ( X , Y ) (cid:54) = ∅ , then there exists a nonnegative equilibrium if and only if β γ > | ( A α · , →− α · , ) , || A , | ; the same is true if positive interfacing is indirect, R + γ = R + γ ( X , Y ) (cid:54) = ∅ , and if | A , | (cid:54) = 0 . (iv) Combined direct negative interfacing . Let R − γ = R − γ ( X , Y ) (cid:54) = ∅ . If positive interfacingis direct, R + γ = R + γ ( X , Y ) (cid:54) = ∅ , then there exists a nonnegative equilibrium for any choice of β ∈ R > and γ ∈ R > ; the same is true if positive interfacing is indirect, R + γ = R + γ ( X , Y ) (cid:54) = ∅ ,and if | A , | (cid:54) = 0 .Proof. Equilibria of the output network, denoted by ( x ∗ , y ∗ , y ∗ ) ∈ R N +2 , satisfy = α · , + A x ∗ + γ y ∗ e i − γ x ∗ j y ∗ e j , for i, j ∈ { , } ,x ∗ = β γ ∈ R > , y ∗ y ∗ = β β , where y ∗ , y ∗ / ∈ { } , (32)where ( α · , + A x ) is the unknown kinetic function from (31), see also Lemma A.1. Statements (iii)and (iv) of the theorem are now proved; statements (i) and (ii) follow analogously.22 ase : R + γ = R + γ ( X , Y ) (cid:54) = ∅ and R − γ = R − γ ( X , Y ) (cid:54) = ∅ . Defining vectors ¯x ∗ ≡ ( x ∗ , x ∗ , . . . , x ∗ N ) (cid:62) , ¯ α · , ≡ ( α , , α , , . . . , α N, ) (cid:62) , ¯ α , · ≡ ( α , , α , , . . . , α ,N ) (cid:62) and ¯ α · , ≡ ( α , , α , , . . . , α N, ) (cid:62) , equa-tion (32) can be written as0 = (cid:18) α , − α , β γ + (cid:104) ¯ α , · , ¯x ∗ (cid:105) (cid:19) y ∗ + γ β β , ¯x ∗ = − ( A , − γ y ∗ diag( e )) − (cid:18) ¯ α · , + ¯ α · , β γ (cid:19) . (33)Assume there exists y ∗ >
0. Since all eigenvalues of A have negative real parts, the same is truefor A , , implying nonpositivity of matrix ( A , − γ y ∗ diag( e )) − ∈ R ( N − × ( N − ≤ [29, Theorem4.3], so that ¯x ∗ ∈ R N − ≥ . Matrix ( A , − γ y ∗ diag( e )) − can be written in a cofactor form( A , − γ y ∗ diag( e )) − | A , | − γ y ∗ | A (1 , , (1 , | = C , C , . . . C N − , C , C , . . . C N − , ... ... . . . ... C ,N − C ,N − . . . C N − ,N − − γ y . . . D , . . . D N − , ... ... . . . ...0 D ,N − . . . D N − ,N − , (34)where C i.j and D i,j are the ( i, j )-cofactors of A , and A (1 , , (1 , , respectively. Substituting (34)into (33), one obtains the quadratic equation0 = γ | A , || A , | (cid:18) β γ − | ( A α · , →− α · , ) , || A , | (cid:19) y − (cid:20) | A || A , | (cid:18) β γ − | A α · , →− α · , || A | (cid:19) − γ γ β β | A (1 , , (1 , || A , | (cid:21) y − γ β β . (35)Since A is cross-nonnegative with a negative spectral abscissa, it follows that | A , | / | A , | > | A | / | A , | < | A (1 , , (1 , | / | A , | < | ( A α · , →− α · , ) , | / | A , | ≥
0, and | A α · , →− α · , | / | A | ≥
0. If( β /γ −| ( A α · , →− α · , ) , | / | A , | ) >
0, it follows from (35) that there exists an equilibrium ( y ∗ , y ∗ ) ∈ R > , consistent with the assumption. On the other hand, if ( β /γ − | ( A α · , →− α · , ) , | / | A , | ) ≤ | A α · , →− α · , | / | A | ≥ | ( A α · , →− α · , ) , | / | A , | , it follows that ( y ∗ , y ∗ ) / ∈ R > . Case : R + γ = R + γ ( X , Y ) (cid:54) = ∅ and R − γ = R − γ ( X , Y ) (cid:54) = ∅ . Equation (32) reads0 = α , − α , β γ + (cid:104) ¯ α , · , ¯x ∗ (cid:105) , ¯x ∗ = − ( A , − γ y ∗ diag( e )) − (cid:18) ¯ α · , + ¯ α · , β γ + γ β β ( y ∗ ) − e (cid:19) , (36)from which, using (34), one obtains0 = γ | A , || A , | (cid:18) β γ − | ( A α · , →− α · , ) , || A , | (cid:19) y − | A || A , | (cid:18) β γ − | A α · , →− α · , || A | (cid:19) y + γ β β | A , || A , | . (37)If | A , | (cid:54) = 0 and ( β /γ − | ( A α · , →− α · , ) , | / | A , | ) >
0, using the fact that | A , | / | A , | <
0, it fol-lows that there exists an equilibrium ( y ∗ , y ∗ ) ∈ R > . If | A , | (cid:54) = 0 and ( β /γ −| ( A α · , →− α · , ) , | / | A , | ) ≤
0, then ( y ∗ , y ∗ ) / ∈ R > . 23 ase : R + γ = R + γ ( X , Y ) (cid:54) = ∅ and R − γ = R − γ ( X , Y ) (cid:54) = ∅ . Equation (32) reads0 = γ ( y ∗ ) + (cid:18) α , − α , β γ + (cid:104) ¯ α , · , ¯x ∗ (cid:105) (cid:19) y − β γ β γ , ¯x ∗ = − A − , (cid:18) ¯ α · , + ¯ α · , β γ (cid:19) ∈ R N − ≥ , (38)from which it follows that ( y ∗ , y ∗ ) ∈ R > for any choice of β ∈ R > and γ ∈ R > . Case : R + γ = R + γ ( X , Y ) (cid:54) = ∅ and R − γ = R − γ ( X , Y ) (cid:54) = ∅ . Equation (32) reads0 = (cid:18) α , − α , β γ + (cid:104) ¯ α , · , ¯x ∗ (cid:105) (cid:19) y ∗ − β γ β γ , ¯x ∗ = − A − , (cid:18) ¯ α · , + ¯ α · , β γ + γ y ∗ e (cid:19) . (39)Assume there exists y ∗ >
0. Then, it follows that x ∗ ∈ R N − ≥ . Substituting (34) with γ = 0into (39), one obtains the quadratic equation0 = γ | A , || A , | ( y ∗ ) − | A || A , | (cid:18) β γ − | A α · , →− α · , || A | (cid:19) y ∗ + β γ β γ . (40)If | A , | (cid:54) = 0, there exists an equilibrium ( y ∗ , y ∗ ) ∈ R > for any choice of β ∈ R > and γ ∈ R > ,consistent with the assumption. Remark . If the x -component of the asymptotically stable equilibrium of the input network R α ( X )is zero, x ∗ = 0, then, like case (iv), cases (i) and (iii) from Theorem E.1 also generically ensurean existence of a nonnegative equilibrium; case (ii) then necessarily leads to negative equilibria.However, stable unimolecular networks with x ∗ = 0 are a negligible subset of the more generalstable unimolecular networks and are, hence, of limited practical relevance.Let us note that Theorem E.1 has an intuitive interpretation. In particular, x ∗ = | A α · , →− α · , | / | A | is the x -component of the equilibrium of the input network, ( x ∗ ) | = | ( A α · , →− α · , ) , | / | A , | isthe x -component of the equilibrium of the restricted input network R α ( X \ X ) with species con-centration x ≡
0. Furthermore, condition | A , | (cid:54) = 0 requires that the species X influences X inthe input network; more precisely, | A , | (cid:54) = 0 requires that there exists at least one directed pathfrom X to X in the digraph induced by the Jacobian matrix of R α ( X ).We now prove that if, for a particular choice of the rate coefficients, the output network fromTheorem E.1 has no nonnegative equilibria, then the network displays deterministic and stochasticblow-ups, i.e. we prove that the output network undergoes deterministic and stochastic NECs(see also Definition A.4 in Section A.2.3). In what follows, the critical values of β /γ at whichnonnegative equilibria cease to exist in Theorem E.1 are called bifurcation points . Theorem E.2. (Negative-equilibrium catastrophe)
Consider a unimolecular input network R α whose RRE s have an asymptotically stable equilibrium, and a controller R β,γ of the form (8) . Then,excluding the bifurcation points, the output network R α,β,γ = R α ∪ R β,γ displays deterministic andstochastic negative-equilibrium catastrophe for all nonnegative initial conditions.Proof. Consider the case with pure direct positive interfacing from Theorem E.1(i); the RREs forthe concentration x = ( x , x , . . . , x N ) ∈ R N ≥ and ( y − y ) ∈ R are given byd x d t = α · , + A x + γ y e , dd t ( y − y ) = γ x − β . (41)24et w ≡ | A | − ( C , , C , , . . . , C N, ) ∈ R N , where C i,j is the ( i, j )-cofactor of matrix A ∈ R N × N .Taking the inner product (cid:104) w , ·(cid:105) in the first equation from (41), and using the fact that (cid:104) A (cid:62) w , x (cid:105) = x , one obtains: dd t (cid:104) w , x (cid:105) = | A α · , → α · , || A | + x + γ | A , || A | y . (42)Using the fact that | A , | / | A | <
0, equations (41) and (42) imply thatdd t (cid:0) −(cid:104) w , x (cid:105) + γ − ( y − y ) (cid:1) = (cid:18) | A α · , →− α · , || A | − β γ (cid:19) − γ | A , || A | y ≥ (cid:18) | A α · , →− α · , || A | − β γ (cid:19) . (43)By Theorem E.1(i), a nonnegative equilibrium does not exist if and only if β γ ≤ | A α · , →− α · , || A | ; ex-cluding the bifurcation point β γ = | A α · , →− α · , || A | , it follows from (43) that the linear combination( −(cid:104) w , x (cid:105) + γ − ( y − y )), and hence an underlying concentration, is a monotonically increasingfunction of time for all nonnegative initial conditions, i.e. the output network displays a determin-istic NEC. Identical argument implies that then the output network displays a stochastic NEC aswell. Deterministic and stochastic NECs for cases (ii) and (iii) from Theorem E.1 are establishedanalogously. References [1] Endy D., 2005. Foundations for engineering biology.
Nature , : 449–453.[2] Del Vecchio, D., Dy, A. J., Qian, Y., 2016. Control theory meets synthetic biology. Journal ofthe Royal Society Interface , ( ): 3–43.[3] Gardner, T. S., Cantor, C. R., Collins, J. J., 2000. Construction of a genetic toggle switch inEscherichia coli. Nature , : 339–342.[4] Elowitz, M. B., Leibler, S., 2000 A synthetic oscillatory network of transcriptional regulators. Nature , : 335–338.[5] Chappell, J., Takahashi, M. K., Lucks, J. B., 2015. Creating small transcription activatingRNAs. Nature chemical biology , ( ): 214–220.[6] Isaacs, F. J., Dwyer, D. J., Ding, C., Pervouchine, D. D., Cantor, C. R., Collins, J. J., 2004.Engineered riboregulators enable post-transcriptional control of gene expression. Nature biotech-nology , ( ): 841–847.[7] Kar S., Baumann W. T., Paul M. R., Tyson J. J., 2009. Exploring the roles of noise in theeukaryotic cell cycle. Proceedings of the National Academy of Sciences
USA, : 6471–6476.[8] Vilar, J. M. G., Kueh, H. Y., Barkai, N., Leibler, S., 2002. Mechanisms of noise-resistance ingenetic oscillators.
PNAS , USA, ( ): 5988–5992.[9] Erban, R., Chapman, J. Stochastic Modelling of Reaction-Diffusion Processes . Cambridge Textsin Applied Mathematics, Cambridge University Press, 2019.2510] Feinberg, M.
Lectures on Chemical Reaction Networks , Delivered at the Mathematics ResearchCenter, U. of Wisconsin, 1979.[11] Gillespie, D. T., 1992. A rigorous derivation of the chemical master equation.
Physica A:Statistical Mechanics and its Applications , (1): 404–425.[12] Drengstig, T., Ueda, H. R., Ruoff, P., 2008. Predicting perfect adaptation motifs in reactionkinetic networks. Journal of Physical Chemistry B , ( ): 16752–16758.[13] Ferrell, J. E., 2016. Perfect and near-perfect adaptation in cell signaling. Cell Systems , ( ):62–67.[14] Chandra, F. A., Buzi, G., Doyle, J. C., 2011. Glycolytic oscillations and limits on robustefficiency. Science , ( ): 187–192.[15] Barkai, N., Leibler, S., 1997. Robustness in simple biochemical networks. Nature , : 913–917.[16] Spiro, P., Parkinson, J., Othmer, H. G., 1997. A model of excitation and adaptation in bacterialchemotaxis. PNAS , USA, : 7263–7268.[17] Yi, T. M., Huang, Y., Simon, M. I., Doyle, J., 2000. Robust perfect adaptation in bacterialchemotaxis through integral feedback control. PNAS , ( ): 4649–4653.[18] Boo, A., Ellis, T., Stan, G. B., 2019. Host-aware synthetic biology. Current Opinion in SystemsBiology , : 66–72.[19] ˚Astr¨om, K. J., and H¨agglund, T. (1995). PID Controllers: Theory, Design, and Tuning .Instrument Society of America, 1995.[20] Plesa, T., Vejchodsk´y, T., and Erban, R., 2016. Chemical Reaction Systems with a HomoclinicBifurcation: An Inverse Problem.
Journal of Mathematical Chemistry , ( ): 1884–1915.[21] Oishi, K., and Klavins, E., 2011. Biomolecular implementation of linear I/O systems. IETSystems Biology , Volume 5, Issue 4: 252–260.[22] Briat, C., Gupta, A., Khammash, M., 2016. Antithetic integral feedback ensures robust perfectadaptation in noisy bimolecular networks.
Cell Systems , ( ): 15–26.[23] Aoki, S.K., Lillacci, G., Gupta, A., Baumschlager, A., Schweingruber, D., and Khammash,M., 2019. A universal biomolecular integral feedback controller for robust perfect adaptation. Nature , : 533–537.[24] Olsman, N., Baetica, A. A., Xiao, F., Leong, Y.P., Doyle, J., and Murray, R., 2019. Hard limitsand performance tradeoffs in a class of antithetic integral feedback networks. Cell Systems , ( ):49–63.[25] Olsman, N., Xiao, F., Doyle, J., 2019. Architectural principles for characterizing the perfor-mance of antithetic integral feedback networks. ISince , : 277–291.[26] Plesa, T., Vejchodsk´y, T., and Erban, R., 2017. Test Models for Statistical Inference: Two-Dimensional Reaction Systems Displaying Limit Cycle Bifurcations and Bistability. StochasticDynamical Systems, Multiscale Modeling, Asymptotics and Numerical Methods for Computa-tional Cellular Biology , 2017. 2627] Gillespie, D.T., 1977. Exact stochastic simulation of coupled chemical reactions.
Journal ofPhysical Chemistry , ( ): 2340–2361.[28] Cox, D., Little, J., and O’Shea, D. Using algebraic geometry . Springer, second edition, 2005.[29] Fiedler, M., Pt´ak, V., 1962. On matrices with non-positive off-diagonal elements and positiveprincipal minors.
Czechoslovak Mathematical Journal , ( ): 382–400.[30] Farina, L., Rinaldi, S. Positive linear systems: Theory and applications . John Wiley & Sons,Inc, 10.1002/9781118033029, 2000.[31] Gillespie, D. T., 1992. A rigorous derivation of the chemical master equation.
Physica A:Statistical Mechanics and its Applications , (1): 404–425.[32] Van Kampen, N. G. Stochastic processes in physics and chemistry . Elsevier, 2007.[33] Kurtz, T. G., 1972. The relationship between stochastic and deterministic models for chemicalreactions.
Journal of Chemical Physics , : 2976–2978.[34] Plesa, T., Stan, G. B., Ouldridge, T. E., and Bae., W., 2021. Quasi-robust control of biochem-ical reaction networks via stochastic morphing. Available as https://arxiv.org/abs/1908.10779.[35] Klonowski, W., 1983. Simplifying principles for chemical and enzyme reaction kinetics. Bio-phys. Chem. (3