Absolutely Robust Controllers for Chemical Reaction Networks
AAbsolutely Robust Controllers for Chemical ReactionNetworks
Jinsu Kim German EncisoJuly 27, 2020
Abstract
In this work, we design a type of controller that consists of adding a specific set of reactionsto an existing mass-action chemical reaction network in order to control a target species. Thisset of reactions is effective for both deterministic and stochastic networks, in the latter casecontrolling the mean as well as the variance of the target species. We employ a type of net-work property called absolute concentration robustness (ACR). We provide applications to thecontrol of a multisite phosphorylation model as well as a receptor-ligand signaling system.For this framework, we use the so-called deficiency zero theorem from chemical reactionnetwork theory as well as multiscaling model reduction methods. We show that the targetspecies has approximately Poisson distribution with the desired mean. We further show thatACR controllers can bring robust perfect adaptation to a target species and are complementaryto a recently introduced antithetic feedback controller used for stochastic chemical reactions.
Keywords absolute concentration robustness, control, reaction networks, Possion distribution,multiscaling, deficiency zero.
In this paper we propose a set of synthetic controllers that can be added to a given chemicalreaction network in order to control the concentration or copy number of a given species of interest.Chemical reaction networks describe a variety of problems in engineering and biology, and therehas recently been a surge in interest for stochastic models of such networks [1–4]. Stochasticeffects are important in order to describe the noise inherent in reactions with low numbers ofmolecules, as is often the case inside individual cells. The techniques proposed in this paper willbe shown to apply both in the deterministic case, where concentrations are described by ordinarydifferential equations, as well as in the stochastic case where the dynamics are described by acontinuous time, discrete space Markov process.The controllers used in our framework are inspired by a property called absolute concentrationrobustness (ACR), which guarantees that a species has the same steady state value regardless ofthe initial conditions. We provide both a theoretical framework and computational simulations inseveral specific biochemical systems to show that an ACR controller can shift all positive steadystate values of a target species towards a desired value. We also show that in stochastic networks1 a r X i v : . [ q - b i o . M N ] J u l atisfying certain topological criteria, the ACR controller can account for the intrinsic noise in thechemical reaction. We approximate the behaviour of the target species using a reduced chemicalreaction model derived through multiscaling analysis. Our stochastic analysis assumes certain con-ditions on the topology of the controlled network that are described using the so-called deficiency and weak reversibility of the system. These two theoretical tools will be combined to calculate thebehaviour of the reduced system, as well as to show that the behaviour of the target species in thereduced system approximates that of the original network. Using computational simulations wealso explore the robust perfect adaptation of the target species in the controlled system, a highlydesirable goal in control theory.When a dynamical system has multistationarity or the dynamics are confined to a lower-dimensional subset by conservation relations among species, the long term behaviour dependson the initial conditions of the system. In general different initial concentrations may lead to dif-ferent long term steady states of the different species. However sometimes the positive steady statevalues of a species of interest are identical independent of the initial conditions. Such a system issaid to possess the ACR property as the steady state of the dynamics is robust to the initial condi-tions. In that case, the species with identical steady state values is called the ACR species. Thiscounter-intuitive dynamical aspect was proposed by Shinar and Feinberg in 2010 [5], where theyfurther provided network topological conditions ensuring that the associated deterministic systemadmits ACR.For a simple example of an ACR system consider the following network, which will be thebasic ACR controller throughout this manuscript, Z + A θ −→ Z, Z µ −→ A. (1)Both reactions produce and consume the same amount of A and Z , hence the total amount of Z + A is conserved. One can think of A and Z as being different forms of the same protein, say active andinactive. Let a ( t ) and z ( t ) be the concentration of the species A and Z , respectively. Assumingthe associated dynamical system is equipped with mass-action kinetics [6], the concentration of Z follows the equation ddt z ( t ) = θa ( t ) z ( t ) − µz ( t ) . (2)At steady state one can set the right hand side equal to zero, and assuming z (cid:54) = 0 one obtains thatthe steady state value for A is µθ . Letting a (0) + z (0) = N be the initial input of the system, theonly positive roots of the right-hand side of (2) are ( a, z ) = ( µθ , N − µθ ) . Hence this system is anACR system and species A is an ACR species. See Figure 1b for a phase plane diagram illustratingthis behaviour.In comparison to the ACR network in (1), we consider the simple reaction A − (cid:42)(cid:41) − B, (3)where the total mass of A and B is also invariant in time. In both systems the dynamics is confinedto one of the black straight lines in Figure 1a and 1b. The positive steady states of this system lieon the intersections between the nullclines (red) and the phase planes (black). Hence for system(3), the steady state values of both A and B vary depending on the initial condition as shown in2igure 1a. On the other hand, the ACR network system (1) is such that all the positive steady statesfor A are identical, as shown in Figure 1b.A similar type of control has been considered by Mustafa Khammash and others in Briat et. al.[7]. In that work, Khammash and colleagues propose an antithetic integral feedback circuit Z + Z η −→ µ −→ Z k −→ Z + X , X (cid:96) θ −→ X (cid:96) + Z , that robustly stabilizes a species of interest, X (cid:96) , in the presence of intrinsic noise. In the controlledsystem the mean of the stochastic dynamics of a target species is stabilized at a pre-specified valuewith a low metabolic cost. A recent follow-up work [8] experimentally implemented the antitheticcontrol circuit in a growth-rate control system in E. coli . We point out that this antithetic controlcircuit satisfies the topological conditions for ACR provided by Shinar and Feinberg in [5], and thusit is a specific example of an ACR system. In particular, the control species Z and Z are ACRspecies if the system admits a positive steady state. There are however a few important differenceswith our work. As we will show, the ACR controllers aim to control the mean, variance and evenhigher moments of a target species by controlling its distribution. The controllers in [7, 8] aredesigned to robustly control the mean of the target species, but the controller might increase itsnoise. To account for the noise in the target species, Briat et. al. show in the follow-up work [9]that for a unimolecular model, an additional negative feedback loop can reduce the noise up to theoriginal variance. Also, stochastic ACR controllers control the target distribution approximatelyassuming that sufficient copies of the control species present, while the control proposed by [7, 8]provides exact control without approximation scheme.A particularly powerful property of an ACR system is that it can endow the ACR property toa given network. When an ACR system is added to a non-ACR model, one of the species in thecombined system could become absolutely robust. For example, suppose that species A is presentin a given deterministic network but that Z is not, and that we add the two reactions (1). Thenthe dynamics of Z will still satisfy z ( t ) = θaz − µz . Moreover, at steady state it must still hold a = µ/θ due to the same analysis as before. Thus species A is now absolutely robust in the newsystem.If the ACR property of the ACR system is inherited to the target species in the controlledsystem, then we call the ACR system an ACR controller . Throughout this paper we also call thenewly introduced species in the ACR controller a control species . In the following sections, weshow that the steady state value of the target species is tunable with the parameters of the ACRcontroller. In the Supplementary Material, we further investigate the local stability of the steadystate in the controlled system.While the ACR controller has the ability to control a target species in deterministic systems,chemical species are often modeled as discrete entities. Stochastic models in biology have becomeincreasingly relevant, as people have noticed that intrinsic noise significantly contributes to thedynamical behaviour [10–15]. The effects of noise are especially large if the abundance of aspecies in the system is low. Many important biochemical models consist of species with low copynumbers inside each individual cell [4]. More details on the modeling of stochastic networks areincluded in the Supplementary Material. In stochastic models we have additional control goalsthan for deterministic models, as it is important to not only control the mean expression level butalso its variance (i.e. noise) and ideally the full probability distribution of the target species. If athermostat makes the temperature of a room oscillate between ◦ F and ◦ F, one might arguethat the room is controlled with mean temperature ◦ F, although its occupants might disagree.3n order to control stochastic systems, we rely on the mathematical theory of deficiency inchemical reaction networks. The deficiency of a reaction network is a non-negative integer thatis determined by the topology of the network regardless of parameter values. Networks with de-ficiency zero and a weak reversibility property have well characterized long term dynamics, un-der both deterministic and stochastic conditions. In the deterministic case, such systems admita unique local asymptotically stable steady state for given total amounts of the species [16–18].For a stochastic system, under the same conditions, each of the species has Poisson distributioncentered around its deterministic steady state [1] (see the Supplementary Material for additionaldetails). These strong properties inspire us to propose a new deficiency based control scheme forstochastic reaction networks, based on recent work expanding ACR to the stochastic case [19, 20].One property observed in some stochastic chemical reaction networks is a so-called extinctionevent . Such an event takes place when some of the species disappear and can never return to thesystem. A stochastically modeled ACR controller can go extinct if a control species, such as Z in the basic ACR controller, is entirely removed from the system. This phenomenon is commonlypresent in ACR networks [21]. One way to minimize this effect is to run the controllers withsufficiently high control species abundance, so that a potential breakdown of the ACR system israre.This high abundance setting is indeed a suitable assumption for the study of stochastic systems[22]. In our stochastic systems, each species can be categorized as either high abundance or lowabundance, compared with the total protein abundance N . We use N as a scaling parameter, andwe carry out a multiscaling procedure to reach a reduced stochastic reaction network. By assumingthat the reduced network has zero deficiency and is weakly reversible, we conclude that the targetspecies has approximately Poisson distribution both in the reduced and the original networks. Inspecial cases, the reduced model can be treated as a hybrid between deterministic and stochasticnetworks [20, 23, 24].We provide multiple examples of control of given biochemical networks using different ACRcontrollers. For example, we use an existing deterministic model of ERK signal transduction fromRubinstein et al [25], and we stabilize the dose response of this system using the basic ACR con-troller. We also study a stochastic receptor-ligand model in which we target the concentration offree receptor, and we show that the concentration of a downstream regulatory protein is also con-trolled as a result. Finally, we study a stochastic dimer-catalyzer model together with an expandedACR controller as an application of the hybrid approach. Simulations using the Gillespie algorithmare provided throughout to illustrate the control implemented by our approach. We also emphasizethat the controlled networks admit a robust perfect adaptation property for both deterministic andstochastic examples. We begin by using an ACR controller for a deterministic system. Consider again the simple trans-lation model between species A and B , A k − (cid:42)(cid:41) − k B. (4)4 Time
Uncontrolled
Total mass=100Total mass=200Total mass=500
A BA + ZZ A κ κ θ µ Original system
ACR controller
AB AZ
A B A +ZZ A2Z a Controlled
A A dc e
Time
Total mass=100Total mass=200Total mass=500
Figure 1: a. and b. Dynamics of the networks A (cid:10) B and the basic ACR controller, respectively.The intersection between each black line and the red line is a steady state for a given total mass. c. The original model in 1a is controlled using the basic ACR controller. d. Time evolution of A inthe original system, k = 2 , k = 4 . e. Time evolution of A in the controlled system via the ACRcontroller, using the above parameters and θ = 1 , µ = 20 .Letting a ( t ) and b ( t ) denote the concentration of species A and B , respectively, the associateddeterministic system with mass-action kinetics is ddt a ( t ) = − k a ( t ) + k b ( t ) , ddt b ( t ) = k a ( t ) − k b ( t ) . (5)We notice that ddt a ( t ) + ddt b ( t ) = 0 which implies that the total mass a ( t ) + b ( t ) is conserved.When a (0) + b (0) = N , the steady state of the system is ( a ∗ , b ∗ ) = ( κ κ + κ N, κ κ + κ N ) by usingthe conservation a ( t ) + b ( t ) = N . Hence the positive steady state concentration of A in theoriginal system (4) varies along with the initial input N . To get the desired steady state value for A , therefore, fine-tuning of the initial condition N is necessary.Adding the basic ACR controller (1) to the original system, we have a new system Z + A θ −→ Z, Z µ −→ A k − (cid:42)(cid:41) − k B. (6)As described in the introduction, using the equation for Z we can deduce that for any positivesteady state it must hold a ∗ = µθ . See Theorem S3.2 in the Supplementary Material for a general-ization of this statement to other networks as well as systems with reaction kinetics different frommass action.As an application, we use the basic ACR network to control a system of an extracellular signalregulated kinase (ERK) activation shown in Rubinstein et al [25]. ERK is a widely studied proteinin signal transduction, and it is activated through phosphorylation at two different sites. The steps5f the dual phosphorylation are regulated by other protein kinases [25]. We denote ERK by S , andconsider the four phosphorylation forms S , S , S and S depending on the phosphorylatedsites. Nonsequential ERK phosphorylation is mediated by mitogen-activated protein kinase MEK,denoted here by E . The variable F denotes a nonspecific phosphatase that mediates ERK dephos-phorylation. The steps of phosphorylation and dephosphorylation are described with the reactionnetwork model in Figure 2a.In the ERK system in Figure 2a, there are three conservation relations. For instance, E tot = E + [ ES ] + [ ES ] + [ ES ] represents the total concentration of kinase, and similarly for totalsubstrate S tot and total phosphatase F tot . It has been shown that the mass action deterministicmodel associated with the ERK model in Figure 2 a has different long-term dynamical behaviourdepending on the system parameters [25–27]. These dynamical behaviours include unique stablestationarity, sustained oscillations, and bistability. We use the parameters in Rubinstein et al.,which are such that the ERK system converges to a unique, stable, and positive steady state [25].One of the most important features of this system is its so-called dose response, which describesactive ERK S as a function of the kinase intput E . However this dose response depends on thetotal amount of phosphatase F tot . We introduce a control using the basic ACR controller in orderto fix S for every given value of total kinase E tot , and therefore to stabilize the dose response.As the plot on the right-hand side of Figure 2b shows, the steady state concentration of protein S is sigmoidal as a function of E tot for fixed F tot . The goal of control with the basic ACR system inFigure 2a is to equalize the positive steady state of the phophatase F for any F tot , and eventuallyto obtain the same sigmoidal curve for the steady state concentration of S , see Figure 2c.Aside from the mathematical model, it is important to think how the basic ACR system in Fig-ure 2a could be implemented experimentally. There are several possible approaches which mightdepend on the individual system. In this case, suppose that the phosphatase F is bifunctional,acting as a phosphatase in its standard form and turning into a kinase Z when it is itself phos-phorylated. Suppose that kinase Z mediates the phosphorylation of protein F as depicted withthe reaction F + Z → Z . Finally, another phosphatase, which is not explicitly modeled in thissystem, eventually dephosphorylates Z into F as described with the reaction Z → F . This set ofassumptions would suffice to implement the control network. Notice that bifunctional enzymes canbe found in the literature, for instance EnvZ in E. coli osmolarity regulation [28]. Notice also thatthe self-mediated phosphorylation can be found in the epidermal growth factor receptor (EGFR)[29, 30]. While the possibility of the practical implementation of the basic ACR controller remainsopen, we focus on theoretic aspects of the controller in this manuscript.Assuming that the dynamics associated with the ERK model and the basic ACR system inFigure 2a follows mass-action kinetics, we use µ = 2 and θ = 1 . This implies that for any input E tot , F tot and S tot , the steady state of F is . The convergence of F to is also theoretically proven,see Section S5.1 in the Supplementary material. Thus in the controlled system, the concentrationof F converges to , unlike the uncontrolled original ERK model which has different steady stateconcentrations of F for different values of F tot , as described in Figure 2b (left) and 2c (left). As aresult, S in the controlled system has identical dose response regardless of the value of F tot (rightplot in Figure 2c). On the other hand, the S dose responses are different in the original system(the right plot in Figure 2c). 6 Uncontrolled
2Z Z F θ µE + S ES ES E + S E + S ES E + S κ κ κ κ κ κ κ κ κ F + S F S F S F + S F + S F S F + S κ κ κ κ κ κ κ κ κ ACR controllerOriginal ERK model ab F S s t e a d y s t a t e Controlled cd e
Time F S s t e a d y s t a t e F tot =20F tot =30F tot =40 F+Z
Uncontrolled Controlled S S F0 10 for 15 ≤ t ≤ Time F tot =20F tot =30F tot =40 E tot F tot =20F tot =30F tot =40 Perturbation by 0 F F tot =20F tot =30F tot =40 E tot Time
Time
Perturbation by 0 F k F TimeTime
Controlled ERK modelUncontrolled ERK model
Time F Figure 2: a. Reaction network of an ERK system originally proposed in [25], together with the basic ACR controller.Parameters used are k = 3 , k = 2 , k = 1 , k = 2 , k = 1 , k = 3 , k = 2 , k = 5 , k = 3 , k = 2 , k =2 , k = 3 , k = 1 , k = 3 , k = 1 , k = 1 , k = 4 , k = 4 for the original ERK system, and θ = 1 , µ = 2 for the ACR controller. b. (Left) Time evolution of F in the original system without the ACR controller. (Right)Dose response of active ERK S as a function of total kinase E tot for the original system. Initial conditions are S = 50 , E = E tot , F = F tot , and zero for all other species. c (Left) Time evolution of F in the controlled system.(Right) Dose response of S as a function of E tot for the controlled system. Initial conditions are S = 50 , Z = 50 , E = E tot , F = F tot , and zero for all other species. d. Robust perfect adaptation of F in the controlled system. We setthe initial conditions as S = 30 , E tot = 30 , F tot = 30 , z (0) = 30 and zero for all other species. e. The dynamics of S in the original (left) or controlled (right) system which is transiently perturbed by reaction −→ F for t ∈ [15 , . .2 Robust Perfect Adaptation One of the main purposes of the control with an ACR system is to create robust perfect adaptation for the target species. A system possesses perfect adaptation if the output of the system returns tothe pre-stimulus level after the value of the input parameter is changed to a different constant level.Furthermore, if the system achieves perfect adaptation independently of the system parameters, itis said to have robust perfect adaptation (RPA) [31–34].Here we show that species F in the controlled ERK model in Figure 2 admits RPA. We persis-tently disturb the parameters by changing them in time such as k as shown in Figure 2d (top). Asexpected, the concentration of phosphatase F converges to different values for each perturbation asdepicted in Figure 2d (middle). However as Figure 2d (bottom) shows, the controlled ERK systemhas robustness to the perturbations for F as its concentration converges to the set-point µθ = 2 regardless of parameter values. This RPA for F basically arises because the steady state of F iscompletely determined by the two reactions in the ACR controller, and it is independent on theparameters k i of the original ERK system. More details are in Section S3 of the supplementarymaterial.In addition to RPA, we also show that the controlled system is robust to a transient change ofthe system structure. We perturb the system by turning on an additional reaction −→ F for time [15 , . The S concentration in the uncontrolled system initially converges, but it immediatelyresponds to the transient in-flow as F is produced during that time interval (Figure 2e (left)). Theconcentration does not return to the previous steady state value after the transient perturbation isturned off. In the controlled system, the phosphorlyated protein S also responds when the tran-sient in-flow is switched on at t = 15 as show in Figure 2e (right). However, it is quickly drivenback to the steady state after the perturbation is switched off at t = 20 . This robustness to the tran-sient structural disturbance stems from the fact that the controlled system admits a single positivesteady state concentration for the target species F regardless of the input states. Hence when thetransient inflow is turned off, F in the controlled system converges to the set-point wherever thecurrent state of the system has been driven by the perturbation. Recall that the basic ACR system (1) consists of the control species Z and the target species A .The fact that Z only directly controls A may impose limitations in some situations. We show inthe following example that an ACR controller with reactions involving other network species canprovide better performance.For example, consider the following reaction network where no conservation relations exist: A κ (cid:37) A + B κ −→ ∅ κ (cid:38) B (7)In this system, two proteins A and B are constantly produced but also degrade each other.Using the parameters κ = 1 , κ = 3 , κ = 5 , the concentration of A decays toward zero as shownin Figure 3d. Despite the addition of the basic ACR system, the concentration of A still decays8o zero as shown in Figure (3)e. See Section S5.2 in the Supplementary Material for additionaldetails about this system, including the existence of positive steady states.We design the expanded controller shown in Figure 3c to include both A and B in the reactions.It can be verified that the mass-action system associated with this controller is ACR, with ACRspecies A . This is because the additional reactions Z + B (cid:10) Z do not change Z , so that theequation for Z is the same as in the base ACR model. Such reactions that have no contribution tothe control species are also used for the antithetic integral controller in [7, 8]. Using θ = 1 , µ =5 , α = 2 , α = 1 , one can see that this controller steers the positive steady state concentration of A to for different initial conditions (Figure 3f). Here the α and α need to be chosen in a certainrange. See Section S4 in the Supplementary Material for more details about control of general2-dimensional systems. a A +B 0 AB κ κ κ A +Z 2ZB +Z Z A θα α µ c A +B 0 AB κ κ κ Original network b A +B 0 AB A +Z 2ZZ A κ κ κ θ µ Original network
ACR controller
Time
Time
Time A d e f A A
Original network
ACR controller
A(0)=5A(0)=10A(0)=15 A(0)=5A(0)=10A(0)=15
Figure 3: a. Original network, using parameters κ = 1 , κ = 3 and κ = 5 . b . The basic ACR system is added tothe original system, with µ = 5 , θ = 1 . c. Expanded ACR controller, with α = 2 and α = 1 . d. The concentrationof A converges to zero in the original system. e. The basic ACR system in 3b fails to control A , as A still convergesto zero for different initial values. f. The concentration of A is driven to the set-point with each initial condition, forthe controlled system in c. When a system contains species with low copy numbers, the intrinsic noise considerably affectsthe system dynamics. Therefore we model the system stochastically using a Markov process.This continuously evolving Markov process defined on a multi-dimensional integer grid has state-dependent transition rates (for more detail see Section S1.1 in the Supplementary material). Inthe context of stochastic control, recent work by Mustafa Khammash and others has proposedcontrolling a target species by adding four reactions. While that framework allows to control themean of the target species, there could be significant variability in its noise. Our ACR approachmakes use of topological properties of the original network to approximate the full distribution ofthe target species.In stochastic networks, if one of the species reaches zero copies, then a subset of the reactionsin the system would be turned off, potentially preventing the species from ever being produced.9uch an extinction event can take place for Z in the basic ACR controller as well as many otherACR systems [21]. In order to avoid this situation, we design the basic ACR system with suffi-ciently high copies of the control species. More generally, we assume that all species are classifiedinto two types: highly abundant species such as control species Z which are of order N for ascaling parameter N , and low abundance species of constant order. We also scale the parametersof the controlled system to make all the reaction propensities of constant order. Under this samescaling, Enciso [19] used the technique of species ‘freezing’ for an ACR system to generate a re-duced network of low abundant species. It was further shown that if the reduced network has zerodeficiency and is weakly reversible, then an ACR species of low order tends to follow a Poissondistribution centered at its ACR value, as time t and the scaling parameter N go to infinity.The work in [19] approximated the distribution of the target species with the help of a reducedstochastic model, which is the limit of the original stochastic network using a multiscaling pro-cedure. Similar types of approaches have been studied using different system scaling, networktopological conditions or state space truncations [2, 3, 20, 35–38]. The multiscaling assumption in[19] is somewhat special in that all reaction propensities have constant order of magnitude up tofinite time.Given a stochastic chemical reaction network, we now add an ACR controller and use the scal-ing procedure described above in order to study the resulting controlled system. To exemplify thiswe consider a model describing the dynamics of a receptor binding to a ligand and generating adownstream response (Figure 4a). Many important biology models involve receptor-ligand inter-actions such as signal transduction, physiological regulation, and gene transcription. In this case aligand L binds to an inactive receptor R on the cell membrane, converting it into an active receptor R . Two active receptors are dimerized, forming the species D which is phosphorylated sequen-tially in three different locations. The triphosphorylated dimer D transmits the signal inside thecell by activating another protein P as shown in Figure 4d. We control the inactive receptor R using the basic ACR system, in order to control the desired amount of active protein P ∗ .Once again, the practical implementation of such a system must depend on the specific recep-tor. We suggest a possible implementation as follows: suppose that a second ligand, called anantagonist, binds to the receptor forming a molecule Z , which prevents the binding of the originalligand (see Figure 4d). Suppose the complex Z facilitates the recruitment of another antagonist toproduce another copy of Z , leading to the reaction Z + R → Z . The reaction Z → R simplyrepresents the natural unbinding of the antagonist from R . Another option could be to think of Z as a misfolded form of R , and of the reaction Z + R → Z as a prion-like effect where amisfolded receptor makes it more likely that a second receptor will misfold. In any of these cases,the introduction of an new molecule into the system (the antagonist or the misfolded protein) leadsto two additional reactions that control the network.10 a b R2R DD D D µZ(0) θ Z(0) κ L(0) κ κ κ κ κ κ κ κ κ Original Receptor-Ligand model R + Z 2ZZ R θ µ ACR controller L ≈ L(0)Z ≈ Z(0)
Reduced model
L + R R2R DD D D D + P κ κ κ κ κ κ κ κ κ κ κ D +P* PP* κ D + P κ κ D + P* PP* d L PPP
PP*
Antagonist ZR D C o p y nu m b e r Ligand, LControl species, Z
Time
Expected number of jumps for LExpected number of jumps for Z
Figure 4: a. Reaction network for the receptor-ligand pathway (green) and the ACR controller (yellow). b. Reducedmodel obtained by freezing L and Z at their initial values, respectively. c. Stochastic time evolution of the copynumbers of L and Z , highlighting the small net change of L and Z in the system by time t = 150 . d. A schematicpicture for the receptor-ligand model and the ACR controller
We let the system start with initial counts L (0) = 1500 , Z (0) = 1000 , R (0) ≤ and theinitial copy numbers of all the other species equal to zero. Hence species L and Z are the highabundance species of order N = 1000 and the other species are of low abundance.As mentioned above, the main idea of the control for this system is to approximate the distri-bution of R by the reduced network in Figure 4b, which we now explain. Parameters are chosenas κ = 0 . × − , κ = 1 . , κ = 1 . , κ = 1 . , κ = 1 . , κ = 1 . , κ = 1 . , κ =1 . , κ = 1 . , κ = 1 . , κ = 1 . , κ = 0 . , θ = 10 − and µ = 5 × − . In order toarrive to this parameter set, parameters κ through κ were randomly chosen in the range [1 , .Parameters κ , θ and µ are associated with reactions involving high abundance species L and Z ,and they were chosen of order N so that the reactions L + R → R, Z + R → Z and Z → R have constant order propensities under mass-action kinetics. Details of the mass-action propensitycomputations are provided in Section S7.1 in the Supplementary Material. Because of the lowpropensities of the reactions relative to N , the expected change of species L and Z by t = 150 aremuch smaller than the copy numbers of L and Z as Figure 4c shows. By neglecting the relativelysmall number of fluctuations for L and Z shown in Figure 4c, we can freeze them at their initial11ounts and obtain a reduced system in Figure 4b. b Uncontrolled Controlled
Uncontrolled Controlled P r o b a b ili t i e s Distribution of R Copy number
Distribution of R
80 85 90 95 100
UnperturbedPerturbed
PerturbedUnperturbed
Uncontrolled Controlled
Distribution of P* Distribution of P*
Copy number Copy number Copy number P r o b a b ili t i e s P r o b a b ili t i e s P r o b a b ili t i e s a Copy number P r o b a b ili t i e s Distribution of R R (0)=20R (0)=50R (0)=10Poisson(5)
60 80 100
Copy number P r o b a b ili t i e s Distribution of P* R (0)=20R (0)=50R (0)=10 Copy number P r o b a b ili t i e s Distribution of R R (0)=10R (0)=20R (0)=50 Copy number P r o b a b ili t i e s Distribution of P* R (0)=10R (0)=20R (0)=50 dc Figure 5:
Gillespie simulations [39] for the distribution of R and activated protein P ∗ . We use the initial values L (0) = 1500 , P (0) = 100 , R (0) ≤ , and the remaining species have zero initial values. For the ACR controller,we set Z (0) = 1000 . a. For the uncontrolled system, distribution at time t = 150 of inactive receptor R (left) andactive protein P ∗ (right). b. (Left) For the controlled system, distribution at time t = 150 of R (left) and P ∗ (right). c and d. Robustness of the system to randomized parameter perturbations and a transient reaction perturbation, setting R (0) = 20 . c. Four distributions of R obtained by simulations of the uncontrolled (left) and controlled (right)receptor-ligand model with randomly perturbed system parameters. d. Unperturbed and perturbed distributions of P ∗ by transiently switched-on reaction −→ R for t ∈ [50 , in the uncontrolled (left) and controlled (right)receptor-ligand model. Recalling that R in the receptor-ligand network in Figure 4a and its associated reduced net-work 4b behave similarly over time, we carry out an analysis of the reduced network. We initiallyignore the reactions D + P → D + P ∗ and P ∗ → P because they only affect the proteins P, P ∗ without an effect on other species. The resulting network is reversible, and it has zero deficiencysince n − (cid:96) − s = 8 − − , where n is the number of complexes, (cid:96) is the number of connected components, and s is the rankof the stoichiometry matrix. The distribution of R in the reduced network converges to a Poissondistribution in the long run [1]. Depending on the number of states of the stochastic system, thelong-term distribution could be a truncated Poisson distribution. However, the reaction in the ACRcontroller are always reduced to the inflow and the outflow of the target species in the reducedsystem. Hence the long-term distribution is always approximated a Poisson distribution. SeeSection S6.3 in the Supplementary Material for a more rigorous statement.12urthermore, the rate of the Poisson distribution associated with R is determined by the steadystate value of the corresponding deterministic system [1]. Since the rate is equal to the mean ofthe Poisson distribution, the mean of R is close to µθ in the long run. Thus species R in thecontrolled system shown in Figure 4a at a sufficiently large finite time t is well approximated bythe Poisson distribution centered at µθ . This is shown in Figure 5b (left) at t = 150 , where for anyinput R the distribution seems almost Poisson( µθ ). Consequently the protein P distribution is alsorobustly stabilized as shown in Figure 5 b (right). On the other hand both mean and variance of R in the original system vary with respect to different inputs (Figure 5a, left), and this causes thedistribution of P ∗ to change accordingly (Figure 5a, right).In an additional analysis, we study the convergence speed of the distribution of the reducedsystem towards a stationary distribution in Section S7.1 of the Supplementary Material. The un-derlying mathematical framework, with an emphasis on the accuracy of the approximation betweenthe controlled network and the reduced system, is further described in our follow-up paper [40].For the receptor-ligand system, the basic ACR module also robustly controls the target speciesto perturbations. We perturb the parameters κ i in Figure 5 using the equation κ (cid:48) i = κ i + r i , wherethe r i are sampled from a uniform distribution on the interval [0 , . As shown in Figure 5c (left),the uncontrolled system generates distinct distributions of R at t = 150 for randomly perturbedparameters in each simulation. On the other hand, the distributions of R at t = 150 generated bythe controlled system with the same parameters closely approximate the Poisson distribution withmean µθ = 5 , as shown in Figure 5c (right).Plots in Figure 5d show how P ∗ robustly behaves with a transient perturbation in the controlledsystem. We perturb the system with a reaction −→ R only for time t ∈ [50 , . Because of thisadditional input, the distribution of P ∗ at t = 150 is shifted to the right for the uncontrolled system(Figure 5d, left). However for the controlled system, Figure 5d (right) shows that its distribution isrobust to the transient perturbation. Recall that in the receptor-ligand system in Section 2.4, the fluctuation of species L and Z in theconcentrations are negligible since the reaction propensities are small compared with their con-centration. However, many classical studies of stochastic systems eliminate this assumption ofsmall reaction propensities, see for instance the classical work by Kurtz [22]. Reaction propensi-ties could also have different orders of magnitude with respect to N . In such cases, the stochasticsystem is modeled under a multiscaling regime, and its behaviour can be studied using a hybriddeterministic-stochastic system [20, 23, 24, 41–44]. In a hybrid system, the counts of some specieschange stochastically while the concentrations of the other species change continuously. We mod-ify the basic controller in order to control such a hybrid system. In this section, using the finitetime stationary distribution approximation in [20], we show that an expanded basic ACR systemcan be used to control a stochastic system under more general scaling.As an example, we provide a dimer-catalyzer model in Figure 6a. In this system the initial copynumber of species X ∗ , X , C, C p and C pp are all of order N = 1000 . Hence using mass actionkinetics all the reactions have order N propensities. Using the framework established by Andersonet al [20], we approximate the original model with the hybrid system in Figure 6d. The stochasticpart of the hybrid system has zero deficiency and is weakly reversible so that the distribution of thetarget species X is Poisson at a finite time t = 5 [1, 20]. The stochastic and deterministic parts are13oupled, as the mean m ( t ) of X at finite time t is determined by the dynamics of the deterministicsystem as depicted in Figure 6e. A flux balance analysis implies that m ( t ) = k x ∗ ( t ) + µz ( t ) k c pp ( t ) + θz ( t ) . InSection S7.2 of the Supplementary Material, we show that m ( t ) converges to the desired value µθ ,as t → ∞ , when the initial concentration of Z is sufficient. The mean m ( t ) actually convergesto µθ quickly as shown in Figure 6e. Therefore unlike the distribution of X in the original systemas shown in Figure 6b, the distribution of X in the controlled system is approximately Poissoncentered at µθ = 5 at time t = 5 for randomly sampled parameters κ i (6c). Time m(t)
Copy number P r obab ili t i e s Controlled distribution of X
UnperturbedPoisson(5)
Original system
ACR controller X ∗ X X + C pp X ∗ + C pp X + X C C p C pp κ κ κ κ κ κ κ κ X 0 X κ c pp (t) + θ z(t) κ x ∗ (t) + µ z(t) κ c(t) κ x (t) Stochastic
Distribution of X m(0.2)m(5) ab cd e
Deterministic X C C p C pp X ∗ + C pp ZX ∗ κ E(X (t)) κ κ κ κ κ κ µ κ E(X(t))
X + Z 2ZX θ µ θ E(X(t)) Z Probability
Uncontrolled distribution of X
Unperturbed
Copy number P r obab ili t i e s Figure 6:
Dimer-catalyzer model with high reaction rates of order N a. The original model and the ACR controller.Parameters are κ = 1 . , κ = 1 . , κ = 1 . , κ = 1 . , κ = 1 . , κ = 1 . , κ = 1 . and κ = 1 . .The parameters are sampled uniformly randomly in [1 , . We use µ = 5 and θ = 1 for the ACR controller. b. Distribution of X at t = 5 in the uncontrolled system in Figure 6a with both the chosen parameters and randomlyperturbed parameters. c. Distribution of X at t = 5 in the controlled system using both the chosen parameters andrandomly perturbed parameters. Red dotted line indicates the Poisson distribution with the rate . d. The controlledsystem is approximated with a hybrid model consisting of stochastic and deterministic parts. e. The mean m ( t ) of thedistribution of X is displayed (blue solid line) along with the distribution of X at times t = 0 . and t = 5 . Discussion
Absolutely robust networks have the property that the steady state value of a target species isindependent of the total mass of the system. In this paper we have provided a class of controllersbased on absolutely robust networks. We define a control species that interacts with the targetspecies, embedding an absolutely robust network into the given network to enforce target speciesrobustness. For deterministically modeled networks, this type of controller not only stabilizesthe target species at the desired value by tuning the parameters of the ACR controller, but alsomakes the species robustly adapted to parameter perturbations and a transiently supplied additionalreactions. We demonstrate control for deterministic system through an ACR system with an ERKmodel. We illustrate some of our results with the so-called base ACR controller, but we show thatother ACR networks can be also be used.We also show that ACR controllers have the ability to control stochastic networks. The needfor control stochastic system is becoming clear in many disciplines of systems and synthetic bi-ology, particularly given the low species counts present in many individual cells. The average ofa species concentration is a deterministic quantity of a stochastic system, thus one might thinkthat a controller used for a typical deterministic system could also implement stochastic control.However for a nonlinear system, studying the dynamics associated with the averages requires non-trivial tools such as moment closure [45]. Even if mean control is valid with a given controller, thesystem may be still out of control if noise is not properly accounted for. Furthermore, because theassociated stochastic system describes molecular counts of each species instead of concentrations,some species might reach a zero state and lead to an extinction event.As a result, for the control of stochastic systems it can be helpful to use advanced mathemat-ical tools such as theoretical analysis of chemical reaction networks. Using an ACR controllerfor stochastic systems here involves two main mathematical tools, multiscaling model reductionand deficiency zero theorems. To avoid a potential breakdown of a controller because of lack ofreactants, we design an ACR controller with high copies of the control species. Using the toolsabove, we show that a species of interest in the controlled system is roughly Poissonian with tun-able mean and variance. Combining the multiscaling model reduction and the zero deficiencycondition, we show that a simple ACR system can control both mean and variance of an inactivereceptor in stochastic receptor-ligand system as the distribution of the inactive receptor roughlyfollows a Poisson distribution centered at the desired value. The controlled stochastic system alsoadmits robust perfect adaptation as does the corresponding deterministic system.We note that the basic ACR controller used throughout this paper has a connection to classicalcontrol theory, as it admits a non-linear integral feedback that is a well-studied characteristic ofrobustly adapted systems [31, 32, 34]. Integral feedback loops arise in many important biologicalphenomena such as bacterial chemotaxis, photoreceptor responses, or MAP kinase activities. Forthe simple mass-action ACR system Z + X θ −→ Z, Z µ −→ X, (8)the concentration of Z satisfies ddt z ( t ) = z ( t )( − µ + θx ( t )) . Dividing by z ( t ) and integrating onboth sides, we obtain log z ( t ) = log z (0) + (cid:90) t ( θx ( s ) − µ ) ds, (9)15hich is a non-linear integral feedback relation. Such types of integral feedback loops appear inmany different biochemical systems [7, 9, 33, 46, 47].One of the major issues on synthetic controllers is the practical implementation of the proposedcontroller. Aoki et al. [8] show that an antithetic controller could be constructed using two controlproteins, σ factor SigW and anti- σ factor RsiW, in an emphE. coli plasmid implementation. For anACR controller, it remains an open question whether its design is practically feasible in vivo or invitro . One key for synthesizing it is the bifunctionality of an enzyme that potentially brings ACRto the system, as it has been observed for other ACR applications [48, 49]. Notice that the controlspecies Z mediates both production and degradation of the target species X in the basic ACRcontroller. We have suggested some ideas for implementing ACR controllers in our examples. Thecontrol species could be obtained by phosphorlyation of a bifunctional target species, by antagonistligand binding, or by a form of protein misfolding.As sufficient network architectural conditions for ACR property have been shown for examplein the regulation of osmolarity in bacteria [5], designing more general ACR controllers could befeasible. Therefore we believe that this new approach introduced in this paper could help controlother biochemical networks in a way that takes into account stochastic effects. Acknowledgements
We would like to thank Eduardo Sontag, Carsten Wiuf, Chuang Xu, Linard Hoessly and JasonDark for key suggestions regarding this work.
Funding
This work is partially supported by NSF grant DMS1763272 and Simons Foundation grant 594598(Qing Nie). This work is also supported by NSF grant DMS1616233.16 upplementary Material
In this section, we provide mathematical models associated with biochemical systems that we usein the main manuscript, starting with the introduction of reaction networks. A biochemical systemcan be described with a reaction network, which consists of constituent species, complexes that arecombinations of species, and reactions between complexes. A triple ( S , C , R ) represents a reactionnetwork where S , C and R are collections of species, complexes and reactions, respectively. Example 4.1.
Consider the following reaction network describing a substrate-enzyme system. S + E (cid:10) SE → E + P, For this reaction network, S = { S, E, SE, P } , C = { S + E, SE, E + P } and R = { S + E → SE, SE → S + E, SE → E + P } . (cid:52) Regarding a reaction network as a directed graph, each connected component is termed a link-age class . A subset Q of complexes in a linkage class is a strongly connected component if andonly if for any two complexes y, y (cid:48) ∈ Q , there exists a path of directed edges connecting from y to y (cid:48) . If every linkage class in a network consists of a single strongly connected component,then the network is weakly reversible . By the definition, in a network ( S , C , R ) , the set of com-plexes C can be decomposed into disjoint linkage classes. Allowing that a single complex can be astrongly connected component, every linkage class is decomposed into disjoint strongly connectedcomponents.For example, for the following network ( S , C , R ) ∅ (cid:10) C, A (cid:10) B → A + B → C (cid:10) B + C, (10)there are two linkage classes {∅ , C } and { A, B, A + B, C, B + C } . Linkage class {∅ , C } consistsof a single strongly connected component. Linkage class { A, B, A + B, C, B + C } has threestrongly connected components { A, B } , { A + B } and { C, B + C } .Each strongly connected component is further classified into two categories. For a stronglyconnected component Q , if there is no path of directed edges connecting from y ∈ Q to y (cid:48) (cid:54)∈ Q ,then Q is a terminal connected component . Otherwise, Q is a non-terminal connected component .A complex contained in a terminal connected component is called a terminal complex , otherwise itis called a non-terminal complex . In (10), strongly connected components {∅ , C } and { C, B + C } are terminal connected components, and the others are non-terminal connected components.We introduce a domain on which the dynamical system associated with a reaction network isdefined. Definition 4.1.
Let ( S , C , R ) be a reaction network. For a x ∈ R d> , we call a set S x = x + span { y (cid:48) − y : y → y (cid:48) ∈ R} ∩ R d> the stoichiometry class. .2 Dynamical systems For a dynamical system of a reaction network, a reaction rate constant κ for each reaction y → y (cid:48) gives a weight on each reaction, and we denote y κ −→ y (cid:48) to incorporate the rate constant. Witha collection of rate constants K , we denote the associated dynamical system for ( S , C , R ) by ( S , C , R , K ) .For mathematical models of reaction networks, we typically assume that the associated sys-tem is spatially well-stirred. In this case the usual mathematical model for a reaction network iseither a system of ordinary differential equation or a continuous-time, discrete-space Markov pro-cess. When each species has high copy number so that intrinsic noise can be averaged out, theconcentration vector x ( t ) of species in a reaction network ( S , C , R ) is typically modeled with adeterministic network system ddt x ( t ) = (cid:88) y → y (cid:48) κ y → y (cid:48) η y ( x ( t ))( y (cid:48) − y ) , (11)where η y : R d> → R > is a rate function associated to a reaction y κ y → y (cid:48) −−−→ y (cid:48) . One of the prevalentchoice of the rate function is mass action kinetics which defines η y ( x ) = x y , where u v = (cid:81) di =1 u v i i for two vectors u, v ∈ R d ≥ .The intrinsic stochasticity of a system is considered when each species in a reaction networksystem has low copy number. For the usual stochastic model, we use a continuous time, discretestate space Markov process X ( t ) ∈ Z d ≥ defined on Z d ≥ = { z ∈ Z d : z i ≥ for each i } . Thetransitions of X are determined by the reaction vectors. Letting h ( t ) be a function such that lim t → h ( t ) t = 0 , the transition probabilities are defined as P ( X ( t + ∆ t ) = z + y (cid:48) − y | X ( t ) = z ) = (cid:88) ¯ y → ¯ y (cid:48) ¯ y (cid:48) − ¯ y = y (cid:48) − y κ ¯ y → ¯ y (cid:48) λ ¯ y → ¯ y (cid:48) ( z )∆ t + h (∆ t ) , (12)where y (cid:48) − y is a reaction vector associated with a reaction y → y (cid:48) , and λ y → y (cid:48) : Z d ≥ → R ≥ is the reaction intensity representing how likely the associated reaction y → y (cid:48) fires.The usual choice of the propensity functions for a stochastic network system ( S , C , R ) is λ y → y (cid:48) ( x ) = x ( y ) , (13)where u ( v ) = d (cid:89) i =1 u i !( u i − v i )! { u i ≥ v i } for u, v ∈ Z d ≥ . This choice of the propensity function is stochastic mass-action kinetics . In both this supplementary material and the main text, we modelboth the deterministic and the stochastic dynamical system under mass action kinetics. Letting K be the set of reaction intensities associated with R , the quadruple ( S , C , R , K ) defines a (eitherdeterministic or stochastic) dynamical system associated with the reaction network ( S , C , R ) .An infinitesimal behavior of the associated X can be described with the infinitesimal generator A [50], A V ( x ) = lim h → E x ( V ( X ( h ))) − V ( x ) h = (cid:88) y → y (cid:48) ∈R λ y → y (cid:48) ( x )( V ( x + y (cid:48) − y ) − V ( x )) , (14)18or a function V : Z d ≥ → R , where E x denotes the expectation of the process whose initial pointis x . The deficiency of a reaction network is a positive integer determined solely by the structure of thenetwork regardless of parameter values. Let ( S , C , R ) be a reaction network with m complexesand (cid:96) linkage classes. Let further s be the rank of the stoichiometric matrix whose i -th column isgiven by i -th reaction in R . The deficiency δ is equal to m − (cid:96) − s. There are a couple of interpretations of the deficiency. First, we can represent the deterministicsystem (11) as ddt x ( t ) = Y A K ψ ( x ( t )) with a stoichiometry coefficient matrix Y , rate constant matrix A K and the rate function φ ( x ) (See[6] for more details). Then the deficiency δ of a network ( S , C , R ) satisfies δ = dim ( Ker ( Y ) ∩ Im ( A K )) . Second, the deficiency roughly stands for redundancy of the network in the following sense. Con-sider the following two networks, ∅ (cid:10) A and ∅ (cid:10) A (cid:10) A. The deficiency of the left reaction network is − − . The deficiency of the right reactionnetwork is − − . This difference stems from the additional reaction A (cid:10) A in theright network. The gain and loss of one A species is already realized with reaction ∅ (cid:10) A . Hencereaction A (cid:10) A is redundant.Zero deficiency combined with weak reversibility of reaction networks implies very strongcharacteristics of the associated system dynamics for both deterministic models and stochasticmodels. Theorem 4.1 (Horn 1972 [16], Feinberg 1972 [17]) . Let ( S , C , R ) be a weakly reversible reactionnetwork with zero deficiency. Then for any choice of rate parameters, the associated deterministicdynamics endowed with the mass-action kinetics admits a unique locally asymptotic stable positivesteady state at each stoichiometry class. The stationary distribution of the associated stochastic process is fully characterized for aweakly reversible network which has zero deficiency.
Theorem 4.2 (Anderson, Craciun and Kurtz 2010 [1]) . Let ( S , C , R ) be a weakly reversible reac-tion network with zero deficiency. Then for any choice of rate parameters, the associated Markovprocess endowed with the stochastic mass-action kinetics admits a stationary distribution, and it is product form of Poissons (or constrained Poissons). That is, for each x ∈ Z d ≥ in the state space,the stationary distribution π satisfies π ( x ) = M d (cid:89) i =1 c x i i x i ! where c = ( c , c , . . . , c d ) is a steady state of the deterministic counterpart and M is the normal-izing constant. For control of a stochastic model, we use Theorem 4.2 to find an approximation of a targetspecies in a controlled system. Details about this procedure is state in Section 9.
In this section we introduce the absolute concentration robustness (ACR) of a reaction network.In order to make use of ACR systems to design a controller, we consider a special class of ACRnetworks, and then we introduce a precise definition of an ACR controller.
Definition 5.1.
Let ˆ x be a solution to the deterministic network system ( ˆ S , ˆ C , ˆ R , ˆ K ) such that ddt ˆ x ( t ) = ˆ f (ˆ x ( t )) . Suppose this system admits a positive steady state. If there exists a species X ∈ ˆ S such that thevalues of X at any positive steady states are all identical, then ( ˆ S , ˆ C , ˆ R , ˆ K ) is called an ACRnetwork system. Furthermore, the species X and the identical positive steady state value of X are called an ACR species and an ACR value, respectively. Especially if the deterministic model isequipped with mass-action kinetics, the system is called a mass-action ACR network system. In some special cases, the ACR property is determined with a single species in a networksystem.
Definition 5.2.
Let ( ˆ S , ˆ C , ˆ R , ˆ K ) be a deterministic network system modeled with ddt ˆ x ( t ) = ˆ f (ˆ x ( t )) If there exists species S i ∈ ˆ S such that { ˆ x : ˆ f i ( x (cid:48) ) = 0 , ˆ x = (ˆ x , . . . , ˆ x d ) ∈ R d> } = { c } for some c > , then the deterministic system is termed a S i -definite ACR system. Remark 5.1. An S i -definite ACR system is an ACR system. For a S i -definite ACR system, anACR species and its ACR value is solely determined by the single equation associated with thespecies S i .A simple mass-action ACR system constructed with only two species is introduced in [5]. Let ( ˆ S , ˆ C , ˆ R , ˆ K ) be mass-action system associated with Z + X θ −→ Z, Z µ −→ X . (15)20ecause any positive roots x ∗ = ( x ∗ , z ∗ ) of the equation ddt z ( t ) = z ( t )( θx ( t ) − µ ) for species Z satisfies x ∗ = µθ , this mass action system is an ACR system. Furthermore since we have { x : z ( θx − µ ) = 0 } = { µθ , x > , z > } , this system is also a Z -definite ACR system by Definition5.2. We termed this system a basic ACR system for X . This ACR system would be mainly usedfor control in the main text.It is shown that there is a broad collection of networks whose associated mass-action system areACR systems. They are characterized using network topological conditions in [5]. In the followingtheorem, e i denotes a vector whose i th entry is one, and the other entries are all zeros. Theorem 5.2 (Shinar and Feinberg 2010 [5]) . Let ( S , C , R ) be a deficiency 1 reaction network.Suppose there are two non-terminal complexes y and ¯ y such that y − ¯ y = ce i for some i ∈ Z > and c (cid:54) = 0 . Then for any set of parameters K , the mass-action deterministic network system ( S , C , R , K ) is a mass-action ACR network system. The controlled deterministic system is basically a union of two deterministic systems; one isa given network system and the other is an ACR system. We formally define the union of twodeterministic network systems. In the definition below, M n,m denote the set of all n × m matricesand I n denotes the n × n identity matrix. Definition 5.3.
Let ( S , C R , K ) and ( ˆ S , ˆ C , ˆ R , ˆ K ) be deterministic network systems modeled with ddt x ( t ) = f ( x ( t )) and ddt ˆ x ( t ) = ˆ f (ˆ x ( t )) , respectively.Let S = { X , . . . , X d , Y , . . . , Y k } and ˆ S = { X , . . . , X d , Z , . . . , Z ˆ k } . Then the union system ofthe deterministic systems ( S , C R , K ) and ( ˆ S , ˆ C , ˆ R , ˆ K ) is a deterministic system such that ddt ¯ x ( t ) = ¯ f (¯ x ( t )) , where ¯ x ( t ) = ( x ( t ) , . . . , x d ( t ) , y ( t ) , . . . , y k ( t ) , z ( t ) , . . . , z ˆ k ( t )) T and ¯ f = Ef + ˆ E ˆ f with E = I d I k ∈ M d + k +ˆ k,d + k , and ˆ E = I d
00 00 I ˆ k ∈ M d + k +ˆ k,d +ˆ k . Now, we define an ACR controller.
Definition 5.4.
Let ( S , C , R , K ) be a deterministic network system and let ( ˆ S , ˆ C , ˆ R , ˆ K ) be an ACRnetwork system such that X ∈ S ∩ ˆ S and ˆ S \ S (cid:54) = ∅ . If the union of the two network systemsis an ACR network system such that X is an ACR species, then ( ˆ S , ˆ C , ˆ R , ˆ K ) is termed an ACRcontroller for ( S , C , R , K ) and the union system is called a controlled system. Furthermore, if ACRcontroller ( ˆ S , ˆ C , ˆ R , ˆ K ) for ( S , C , R , K ) is a mass-action system, then it is termed a mass-actionACR controller for ( S , C , R , K ) . Steady states and stability using an ACR controller
In this section, we show that for any deterministic system modeled with general kinetics, an ACRcontroller endows ACR to the given system and drives the long-term behavior of a target speciestowards the desired value. For the basic ACR controller, the existence of the steady states willnow be verified together with their stability. In this manuscript, every chemical reaction networkis modeled under mass action kinetics. It is notable, however, that Lemma 6.1, Lemma 6.2, andTheorem 6.3 hold for a class of non-mass action networks, for instance under the generalizedkinetics defined in [6] (Definition 2.2).
Lemma 6.1.
Let ( S , C , R , K ) be a deterministic network system such that X ∈ S . Let ( ˆ S , ˆ C , ˆ R , ˆ K ) be a Z -definite ACR system such that X ∈ ˆ S and Z (cid:54)∈ S . If the union system of ( S , C , R , K ) and ( ˆ S , ˆ C , ˆ R , ˆ K ) admits a positive steady state, then it is an ACR system, and X is an ACR species.Proof. Since Z (cid:54)∈ S , the equation ddt z ( t ) = 0 for Z in the union system is same as the equation for Z in ACR system ( ˆ S, ˆ C , ˆ R , ˆ K ) . At the expense of abusing the notation, we let ddt z ( t ) = f z ( x ( t )) and ddt ¯ z ( t ) = f z (¯ x ( t )) be the equations for Z in ( ˆ S , ˆ C , ˆ R , ˆ K ) and the union system, respectively.By definition of the Z -definite ACR system, if f z ( x ) = 0 then there exists a positive real number c such that x = c . Hence, for each positive steady state ¯ x ∗ in the union system, ¯ x ∗ = c and therefore X is an ACR species in the union system.For a given network system ( S , C , R ) , suppose Z (cid:54)∈ S and X ∈ S . Since the basic ACRcontroller (15) for X is a Z -definite ACR system, it is a mass action ACR controller for ( S , C , R ) .We call this basic ACR system the basic ACR controller interchangeably.Lemma 6.1 guarantees that the values of X must be c at any positive steady states as longas a positive steady state exists in the union system. The following lemma provides a sufficientcondition of a given network system ( S , C , R , K ) for existence of a positive steady state in theunion system of ( S , C , R , K ) and the basic ACR controller. We show that if the desired controlvalue is within the range of the observations, then one can control for that value. Lemma 6.2.
Let ( S , C , R , K ) be a deterministic network system modeled with ddt x ( t ) = f ( x ( t )) . Let
P S = { x : x = ( x , x , . . . , x d ) ∈ R d> , f ( x ) = 0 } . Suppose there exists an x ∗ ∈ P S such that x ∗ = µθ , then the basic ACR network system (15) is a mass-action ACR controller for ( S , C , R , K ) , and the controlled system admits a positive steady state. In the following proof, the concatenation w = ( u, v ) for u ∈ R d and b ∈ R denotes a vectorin R d +1 such that w i = u i for i = 1 , , . . . , d and w d +1 = v . Proof.
Let S = { X , . . . , X d } . Let ¯ x = ( x, z ) for each x ∈ R d ≥ and z ∈ R ≥ be a solution to theunion system of ( S , C , R , K ) and the basic ACR network system (15). Then ¯ x satisfies ddt ¯ x ( t ) = ¯ f (¯ x ( t )) , (16)22or some ¯ f . By the construction of the union system, we have ¯ f i (¯ x ) = f i ( x ) − z ( θx − µ ) if i = 1 ,z ( θx − µ ) if i = d + 1 ,f i ( x ) , otherwise . (17)Let x ∗ be a positive steady state of ( S , C , R , K ) such that x = µθ . For any positive value z ∗ , wehave ¯ f (¯ x ∗ ) = 0 where ¯ x ∗ = ( x ∗ , z ∗ ) .For a given choice of θ and µ , it is a necessary condition that there exists x ∗ ∈ P S such that x ∗ = µθ for a given reaction system because unless Z = 0 , the control species Z is only stabilizedwhen X = µθ . We demonstrate this with the following example. Example 6.1.
Consider this system suggested by one of the reviewers of this manuscript: A −→ B (cid:45) (cid:46) Note that this system admits a unique equilibrium (1 , . Hence for the choice of θ = 1 and µ = 2 ,there is no x ∗ in P S such that x ∗ = µθ . Notably the union of the system and the basic ACR system A + Z θ −→ Z and Z µ −→ A does not admit a positive steady because no positive values a ∗ and z ∗ satisfy both − a ∗ + 1 + z ∗ (2 − a ∗ ) , and z ∗ ( a ∗ − . The convergence to positive steady states in general controlled systems with a ACR controlleris more delicate problem since the actual network structure and parameters need probably to bespecified. However, if linear stability condition is held for a given system as well as the conditionsin Lemma 6.2 with some additional conditions, then the controlled system with the basic ACRsystem (15) admits linear stability. Linear stability of a steady state holds if each eigenvaluesof the Jacobian of a dynamical system at the steady state has a strictly negative real part. Thisimplies the dynamical system asymptotically converges to the steady state if its initial state wasclose enough to the steady state.Remark that in case a given system has no conservation relation, the dynamics is not confinedinto a lower dimensional stoichiometry class. Hence if we assume linear stability of the givensystem at a positive steady state x ∗ , all eigenvalues of the Jacobian at the steady state have strictlynegative real parts. Hence we can maintain the linear stability after we add a ACR controller ifthe parameters of the ACR controller are small enough. This is by the fact that the roots of thecharacteristic polynomial are continuous with respect to the coefficients, hence the eigenvalues ofthe Jacobian of the controlled system still have strictly negative real parts.Hence we investigate the stability of the controlled system when a given system ( S , C , R , K ) admits conservation relations. Let x ( t ) = ( x ( t ) , . . . , x d ( t )) be the deterministic model associated23ith ( S , C , R , K ) such as (11) in R d> . Suppose that u , . . . , u k are positive vectors such that u i · ddt x ( t ) = 0 for all t and for each i , where · means the canonical inner product between twofinite dimensional euclidean vectors. This implies that for a fixed initial state x (0) , there exist M i ’ssuch that u i · x ( t ) = M i for all t. Without loss of generality, we suppose u i are linear independent. Then in the following way,we can reduce the system onto a lower dimension system that admits no conservative relations.First note that since we assume the linear independence of u i ’s, we have k ≤ d . Hence usingGaussian elimination and by rearranging the coordinate of x , we have (cid:2) U | I (cid:3) x ( t ) x ( t ) ... x ¯ d ( t ) x ¯ d +1 ( t ) x d ( t ) = M M ... M k , (18)where ¯ d = d − k , the matrix I is the k dimensional identity matrix, U is some ¯ d × k matrix and M i ’s are some constants. Hence we have x ¯ d +1 ( t ) = M − ¯ d (cid:88) i =1 u i x i ( t ) ,x ¯ d +2 ( t ) = M − ¯ d (cid:88) i =1 u i x i ( t ) , ... x d ( t ) = M k − ¯ d (cid:88) i =1 u i x i ( t ) . (19)This implies the variables x ¯ d +1 , . . . , x d are completely determined by relations (19). Then we havethe following reduced system, ddt x i ( t ) = g i ( x ( t ) , x ( t ) , . . . , x ¯ d ( t )) where ,g i ( x , x , . . . , x ¯ d ) = f i (cid:32) x , x , . . . , x ¯ d , M − ¯ d (cid:88) i =1 u i x i ( t ) , . . . , M k − ¯ d (cid:88) i =1 u ki x i ( t ) (cid:33) . (20)for i = 1 , , . . . , ¯ d . Note that this reduced system is specified with the choice of initial state x (0) asthe initial condition determines the conservative quantity M i ’s. Note further that since the steadystate values of x ¯ d + i for i = 1 , . . . , k are completely determined by the steady state values x i for i = 1 , , . . . , ¯ d , the stability of x ( t ) is also determined by the reduced system ( x ( t ) , . . . , x ¯ d ( t )) .The linear stability of the reduced system is investigated with the eigenvalues of Jacobian. We24enote J ( x ∗ ) be the Jacobian of this reduced system at x ∗ , where we abuse the notation since x ∗ isa state in the original system but the Jacobian is for the reduced system.Now we suppose that species X is the control target with the basic ACR controller X + Z θ −→ Z, Z µ −→ X . (21)Suppose that S is involved in at least one conservation relation. Without loss of generality, suppose u = 1 . Then we have new conservation relations in the union system of ( S , C , R , K ) and the basicACR system. We let ¯ M i = (cid:40) u i · x (0) + u i z (0) = M i + z (0) if u i (cid:54) = 0 ,u i · x (0) + u i z (0) = M i , otherwise . (22)Then the new conservative relations are represented as (cid:2) U | u | I (cid:3) x ( t ) x ( t ) ... x ¯ d ( t ) z ( t ) x ¯ d +1 ( t ) x d ( t ) = ¯ M ¯ M ... ¯ M k , (23)where v is the first column vector of U , and U and I are the same matrices as (18). The definitionof u basically means that the control species Z is involved in the same conservation relation as X in the original system ( S , C , R , K ) . Hence the dynamics ¯ x ( t ) = (¯ x ( t ) , . . . , ¯ x d ( t ) , z ( t )) associatedwith the union system can also be reduced to ddt ¯ x i ( t ) = h i ( x ( t ) , x ( t ) , . . . , x ¯ d ( t ) , z ( t )) . (24)where, h i ( x , x , . . . , x ¯ d , z ) = f x , x , . . . , x ¯ d , ¯ M − j = ¯ d (cid:88) i u j x j ( t ) − z ( t ) , . . . , ¯ M k − ¯ d (cid:88) j =1 u kj x j ( t ) − z ( t ) − z ( θx − µ ) , if i = 1 ,z ( θx − µ ) if i = ¯ d + 1 ,f i (cid:32) x , x , . . . , x ¯ d , ¯ M − i = ¯ d (cid:88) i u i x i ( t ) − v z ( t ) , . . . , ¯ M k − ¯ d (cid:88) j =1 u ij x j ( t ) − u i z ( t ) (cid:33) , otherwiseNote that by (22), we have ∂ i h j ( x ∗ , . . . , x ∗ ¯ d ) = ∂ i g j ( x ∗ , . . . , x ∗ ¯ d ) for i, j ∈ { , , . . . , ¯ d } . Then25he jacobian ¯ J ( x ∗ , z ∗ ) for this system at ( x ∗ , z ∗ ) where x ∗ = µθ is ¯ J ( x ∗ , z ∗ ) = ∂ h ( x ∗ , z ∗ ) − θz ∗ ∂ h ( x ∗ , z ∗ ) · · · ∂ ¯ d h ( x ∗ , z ∗ ) ∂ z h ( x ∗ , z ∗ ) ∂ h ( x ∗ , z ∗ ) ∂ h ( x ∗ , z ∗ ) · · · ∂ ¯ d h ( x ∗ , z ∗ ) ∂ z h ( x ∗ , z ∗ ) ... ∂ h ¯ d ( x ∗ , z ∗ ) ∂ h ¯ d ( x ∗ , z ∗ ) · · · ∂ ¯ d h ¯ d ( x ∗ , z ∗ ) ∂ z h ¯ d ( x ∗ , z ∗ ) θz ∗ · · · = ∂ g ( x ∗ ) − θz ∗ ∂ g ( x ∗ ) · · · ∂ ¯ d g ( x ∗ ) ∂ z h ( x ∗ , z ∗ ) ∂ g ( x ∗ ) ∂ g ( x ∗ ) · · · ∂ ¯ d g ( x ∗ ) ∂ z h ( x ∗ , z ∗ ) ... ∂ g ¯ d ( x ∗ ) ∂ g ¯ d ( x ∗ ) · · · ∂ ¯ d g ¯ d ( x ∗ ) ∂ z h ¯ d ( x ∗ , z ∗ ) θz ∗ · · · (25)As the stability of x ( t ) is determined by the stability of its reduced system, we consider thelinear stability of the reduced system (24) of ¯ x ( t ) to study the stability of the union system of ( S , C , R , K ) and the basic ACR controller (21). For the linear stability, we show that the character-istic function for ¯ J ( x ∗ , z ∗ ) will be the characteristic function of J ( x ∗ ) with some perturbation. Toshow this, we use the conventional notation for deterministic | A | for a square matrix A . I denotesan identity matrix, and it could denote a different dimensional identity matrix according to thecontent. We will also use the column/row expansion of deterministic. We further also use the rowdecomposition for determinant. The row decomposition of the determinant means that when A isa square matrix such that the first row A is equal to A (cid:48) + A (cid:48)(cid:48) with some row vectors A (cid:48) and A (cid:48)(cid:48) ,we have | A | = | A (cid:48) | + | A (cid:48)(cid:48) | where A (cid:48) and A (cid:48)(cid:48) are square matrices whose first row is replaced with A (cid:48) and A (cid:48)(cid:48) , respectively. 26hen | λI − ¯ J ( x ∗ , z ∗ ) | = (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ − ∂ g ( x ∗ ) + θz ∗ − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) − ∂ z h ( x ∗ , z ∗ ) − ∂ g ( x ∗ ) λ − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) − ∂ z h ( x ∗ , z ∗ ) ... − ∂ g ¯ d ( x ∗ ) − ∂ g ¯ d ( x ∗ ) · · · λ − ∂ ¯ d g ¯ d ( x ∗ ) − ∂ z h ¯ d ( x ∗ , z ∗ ) − θz ∗ · · · λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ − ∂ g ( x ∗ ) + θz ∗ − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) − ∂ g ( x ∗ ) λ − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) ... − ∂ g ¯ d ( x ∗ ) − ∂ g ¯ d ( x ∗ ) · · · λ − ∂ ¯ d g ¯ d ( x ∗ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ( − ¯ d +2 θz ∗ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) − ∂ z h ( x ∗ , z ∗ ) λ − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) − ∂ z h ( x ∗ , z ∗ ) ... − ∂ g ¯ d ( x ∗ ) · · · λ − ∂ ¯ d g ¯ d ( x ∗ ) − ∂ z h ¯ d ( x ∗ , z ∗ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) λ − ∂ g ( x ∗ ) − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) − ∂ g ( x ∗ ) λ − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) ... − ∂ g ¯ d ( x ∗ ) − ∂ g ¯ d ( x ∗ ) · · · λ − ∂ ¯ d g ¯ d ( x ∗ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) + θz ∗ λ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) · · · − ∂ g ( x ∗ ) λ − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) ... − ∂ g ¯ d ( x ∗ ) − ∂ g ¯ d ( x ∗ ) · · · λ − ∂ ¯ d g ¯ d ( x ∗ ) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ( − ¯ d +2 θz ∗ (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) − ∂ z h ( x ∗ , z ∗ ) λ − ∂ g ( x ∗ ) · · · − ∂ ¯ d g ( x ∗ ) − ∂ z h ( x ∗ , z ∗ ) ... − ∂ g ¯ d ( x ∗ ) · · · λ − ∂ ¯ d g ¯ d ( x ∗ ) − ∂ z h ¯ d ( x ∗ , z ∗ ) . (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (26)Notice that the first term in (26) is equal to λ | λI − J ( x ∗ ) | . We denote by λG ( λ ) , θz ∗ λH ( λ ) and ( − ¯ d +2 θz ∗ H ( λ ) the first, the second and the third term in (26), respectively. Hence we have | λI − ¯ J ( x ∗ , z ∗ ) | = λG ( λ ) + θz ∗ λH ( λ ) − ( − ¯ d +2 θz ∗ H ( λ ) . (27)Now, using the same notations above, we state a theorem related to the stability of ( x ∗ , z ∗ ) of theunion system of ( S , C , R , K ) and the basic ACR system (21). Theorem 6.3.
Suppose the conditions in Lemma 6.2 hold. Suppose further1. the associated system for ( S , C , R , K ) admits conservative relations such as (18) and u (cid:54) = 0 ,2. all the eigenvalues of J ( x ∗ ) have strictly negative real parts, and . H (0) > if ¯ d is odd and H (0) < if ¯ d is even.Then for sufficiently small θ and µ and for sufficiently large z (0) , all the eigenvalues of ¯ J ( x ∗ , z ∗ ) have also strictly negative real parts. That is the positive steady state ( x ∗ , z ∗ ) is linear stable inthe union system of ( S , C , R , K ) and the basic ACR system (21) . Remark 6.4.
Note that the entries − ∂ z h i ( x ∗ , . . . , x ∗ ¯ d , z ∗ ) can be calculated by using solely f i andthe conservative relations in the original system. For i = 2 , , . . . , ¯ d , by the chain rule ∂ z h i ( x , . . . , x ¯ d , z )= − k (cid:88) j =1 u j ∂f i ∂x ¯ d + j (cid:32) x , x , . . . , x ¯ d , M − ¯ d (cid:88) j =1 x j − z, . . . , M k − ¯ d (cid:88) j =1 u kj x j ( t ) − u kj z ( t ) (cid:33) . Similarly for i = 1 , ∂ z h ( x , . . . , x ¯ d , z )= − k (cid:88) i =1 u j ∂f i ∂x ¯ d + j (cid:32) x , x , . . . , x ¯ d , M − ¯ d (cid:88) j =1 x j − z, . . . , M k − ¯ d (cid:88) j =1 u kj x j − u kj z (cid:33) − ( θx − µ ) . Hence, especially for ( x ∗ , z ∗ ) such that x ∗ = µθ , we have − ∂ z h ( x ∗ , z ∗ ) = − k (cid:88) j =1 u j ∂f i ∂x ¯ d + j (cid:32) x ∗ , x ∗ , . . . , x ∗ ¯ d , M − ¯ d (cid:88) j =1 x ∗ j − z ∗ , . . . , M k − ¯ d (cid:88) j =1 u kj x ∗ j − u kj z ∗ (cid:33) . Proof.
First of all, we scale µ = (cid:15) ¯ µ, θ = (cid:15) ¯ θ and z (0) = M(cid:15) for some R . Note that u (cid:54) = 0 byhypothesis 1 in the statement. Then by (22) and (23), we have z ∗ = ¯ M − x ∗ − ¯ d (cid:88) j =2 u j u x ∗ j = ¯ M − µθ − ¯ d (cid:88) j =2 u j u x ∗ j = z (0) + ¯ d (cid:88) j =1 u j u x j (0) − µθ − ¯ d (cid:88) j =2 u j u x ∗ j , for each i such that u i (cid:54) = 0 . Thus for each (cid:15) , we have θz ∗ = (cid:15) ¯ θ (cid:32) R + (cid:15) ¯ d (cid:88) j =1 u ij v x j (0) − (cid:15) µθ − (cid:15) ¯ d (cid:88) j =2 u ij v x ∗ j (cid:33) > by taking sufficiently large R = R ( (cid:15) ) . We denote c ( (cid:15) ) = θz ∗ , then lim (cid:15) → c ( (cid:15) ) = 0 .28y the hypothesis, all roots of G ( λ ) have strictly negative real parts. Let λ , λ , . . . , λ ¯ d bethe roots of λG ( λ ) where λ = 0 and λ i are non-zero roots with strictly negative real parts. Letfurther denote λ ( (cid:15) ) , λ ( (cid:15) ) , . . . , λ ¯ d ( (cid:15) ) the roots of | λI − ¯ J ( x ∗ , z ∗ ) | = λG ( λ ) + c ( (cid:15) ) λH ( λ ) − ( − ¯ d +2 c ( (cid:15) ) H ( λ ) . By the continuity of roots of a polynomial with respect to the coefficients, wehave lim (cid:15) → | λ i ( (cid:15) ) − λ i | = 0 for i = 1 , , . . . , ¯ d . Hence with small enough (cid:15) , we could make thereal parts of λ i ( (cid:15) ) is still negative for each i = 1 , . . . , ¯ d .We turn to show that λ ( (cid:15) ) has also a strictly negative real part. Note that | λI − ¯ J ( x ∗ , z ∗ ) | (cid:12)(cid:12)(cid:12) λ =0 =( λ − λ ( (cid:15) ))( λ − λ ( (cid:15) )) · · · ( λ − λ ¯ d ( (cid:15) )) . Hence ( − ¯ d +1 λ ( (cid:15) ) λ ( (cid:15) ) · · · λ ¯ d ( (cid:15) ) = − ( − ¯ d +2 c ( (cid:15) ) H (0) . (28)Suppose λ ( (cid:15) ) is a complex number with non-zero imaginary part. Then it must be a conjugateof λ i ( (cid:15) ) for some i . Since we choose (cid:15) small enough so that the real part of each λ i ( (cid:15) ) is strictlynegative, λ ( (cid:15) ) has a negative real part. Now we suppose that λ ( (cid:15) ) is a real number. In thiscase, because of the negative real parts of λ i ( (cid:15) ) , the product (cid:81) ¯ di =1 λ i ( (cid:15) ) is negative if ¯ d is odd andis positive otherwise. Thus by the hypothesis 2 and (28), we have λ ( (cid:15) ) < . Thus the resultfollows.It is notable that the condition of sufficient amount of initial Z in Theorem 6.3 is known to beoften satisfied in practical applications. In this section, we introduce an ACR system that controls both a target species and other speciesin a given network system. By using this type of expanded ACR system, we show that a 2-dimensional reaction system, which admits no positive steady states, can be controlled as weshowed in Section 2.3 of the main text. For a two dimensional system with species A and B ,let A be the target species we desire to control. Then we define a mass action ACR system (29)that controls the other species B as well as the target species A , A + Z θ −→ Z, Z µ −→ A, B + Z α − (cid:42)(cid:41) − α Z. (29) Lemma 7.1.
Let ( S , C , R , K ) be a network system such that S = { A, B } . Let x ( t ) = ( a ( t ) , b ( t )) be the associated deterministic system such that ddt a ( t ) = f ( a ( t ) , b ( t )) , and ddt b ( t ) = f ( a ( t ) , b ( t )) . Suppose that there is no positive steady state for ( S , C , R , K ) . Suppose further that for any positiveconstant c , there exists d > such that f ( c, d ) = 0 . Then the union system of ( S , C , R , K ) and (29) admits a positive steady state ( a ∗ , b ∗ , z ∗ ) such that a ∗ = µθ for any µ and θ provided α b ∗ > α when f ( µθ , b ∗ ) > and α b ∗ < α when f ( µθ , b ∗ ) < .Proof. First, the reactions B + Z (cid:10) Z do not change the concentration of Z . Hence A is stillan ACR species as we showed around (15). Let ¯ x = (¯ a, ¯ b, z ) be the solution to the deterministic29ystem associated with the union system of ( S , C , R , K ) and (29). Then we have ddt ¯ a ( t ) = ¯ f (¯ a, ¯ b, z ) = f (¯ a, ¯ b ) − z ( θ ¯ a − µ ) ,ddt ¯ b ( t ) = ¯ f (¯ a, ¯ b, z ) = f (¯ a, ¯ b ) − z ( α ¯ b − α ) , and ddt z ( t ) = z ( θ ¯ a − µ ) . By the hypothesis, there exists a positive constant x ∗ such that ¯ f ( µθ , x ∗ , z ) = f ( µθ , x ∗ ) = 0 forany z . Then if we choose z ∗ = f ( µθ , x ∗ ) α x ∗ − α , we have ¯ f ( µθ , x ∗ , z ∗ ) = 0 . Since ( S , C , R , K ) does notadmits a positive steady state, f ( µθ , x ∗ ) (cid:54) = 0 . Thus z ∗ (cid:54) = 0 .It follows that z ∗ > because we assumed that α x ∗ > α if f ( µθ , x ∗ ) > and α x ∗ < α if f ( µθ , x ∗ ) < . Therefore ( µθ , x ∗ , z ∗ ) is a positive steady state of the union system. Remark 7.2.
In the proof above, it was shown that b ∗ is independent of α and α . Hence wecan evaluate the value of f ( µθ , b ∗ ) and then we set the parameters α and α in the controller (29)according to the sign of f ( µθ , b ∗ ) . Remark 7.3.
For arbitrary θ > and µ > , if there exists a positive steady state for the controlledsystem ¯ a ( t ) , ¯ b ( t ) and z ( t ) , then it must hold ¯ f ( µθ , b ∗ , z ∗ ) = f ( µθ , b ∗ ) = 0 for some b ∗ . Hence forthe controllablility of the basic ACR system, it is a necessary condition of a 2-dimensional systemthat for any c > there exists d > such that f ( c, d ) = 0 . Remark 7.4.
Stability analysis of the union system with an extended ACR controller is a openproblem.For a higher-dimensional system, we can use a similar idea to build an extended ACR circuitwith multiple control species.
Example 7.1.
Consider the following -dimensional mass-action system admitting no positivesteady states, where protein B and protein C are dimerized to a protein A . B + C −→ A, B ←− −→ C ↑ A Then, we can consider an extended ACR system with two control species Z and Z such that A + Z θ −→ Z , Z µ −→ A,A + Z θ −→ Z , Z µ −→ A,B + Z α − (cid:42)(cid:41) − α Z , C + Z α − (cid:42)(cid:41) − α Z . Let µθ = 10 be the desired set point for A in the union of the given -dimensional system and theextended ACR system. Then species A , B and C are stabilized at (10 , b ∗ , c ∗ ) such that b ∗ c ∗ = 10 ,as long as the parameters α , α , α and α satisfy α > b ∗ α and α > b ∗ α .
10 2051015
A(0)=5A(0)=15A(0)=10
Time
Figure 7: Time evolution of the concentration of A in the union system. In this section, we analyze the stability of the controlled ERK system shown in Figure 2 a. Weshow that both the conditions in Lemma 6.2 and conditions in Theorem 6.3 hold for a positivesteady state ( x ∗ , z ∗ ) in the controlled ERK system. Let ( S , C , R , K ) be the deterministic systemassociated with the ERK network in Figure 2 a using the parameters given in the main text. Letalso ( ¯ S , ¯ C , ¯ R , ¯ K ) be the union system of ( S , C , R , K ) and the ACR controller in 2 a with θ = 1 and µ = 2 . We use a Matlab simulation to obtain a positive steady state, as well as the relevantJacobian, eigenvalues and determinant.Let x ( t ) = ( x ( t ) , x ( t ) , . . . , x ( t )) and ¯ x ( t ) = (¯ x ( t ) , ¯ x ( t ) , . . . , ¯ x ( t ) , z ( t )) represent theconcentrations of species in the systems ( S , C , R , K ) and ( ¯ S , ¯ C , ¯ R , ¯ K ) , respectively. We arrange x , x , . . . , x so that they represent F, E, S , S , S , ES , ES , F S , F S , S , ES and F S , respectively. We also let ¯ x i represent the same concentration as x i , and we let z representthe concentration of Z .To show the condition in Lemma 6.2, we show that the system ( S , C , R , K ) admits a positivesteady state at x ∗ such that x ∗ = µθ = 2 . We verify the existence of the positive steady state usingthe simulation shown in the figure below for the ERK system with F tot = 31 , E tot = 37 , S tot = 100 .31
10 20 30 40 50
Time F F tot =31, E tot =37 and S =100 The positive steady state x ∗ = (2 . , . , . , . , . , . , . , . , . , . , . , . . Werounded off the values to one decimal place.We also notice there are three linear independent conservative relations in ( S , C , R , K ) , x ( t ) = F tot − x ( t ) − x ( t ) − x ( t ) ,x ( t ) = E tot − x ( t ) − x ( t ) − x ( t ) ,x ( t ) = S tot − (cid:88) i =3 x i ( t ) . (30)The target species F is involved in the first conservative relation in (30). Hence condition 1 inTheorem 6.3 holds.For the second condition in Theorem 6.3, we have the following Jacobian J ( x ∗ ) of the reducedsystem obtained by using the conservation laws (30) for ( S , C , R , K ) as (19). J ( x ∗ ) = − . − . . − . . . . − . − . − . − . − . − . . . . − . − . − . . . − . − . − . − . − . .
00 56 . . − . . . . − . . . − . . . − . − . (31)The eigenvalues of J ( x ∗ ) are − . , − . , − . , − . , − . , − . , − . , − . , − . .For the third condition in Theorem 6.3, we note that d = 12 , k = 3 and hence ¯ d = 9 . Thus, if H (0) > , then the condition holds. Note that we have three conservation relations for the system ( ¯ S , ¯ C , ¯ R , ¯ K ) x ( t ) = F tot − x ( t ) − x ( t ) − x ( t ) − z ( t ) ,x ( t ) = E tot − x ( t ) − x ( t ) − x ( t ) ,x ( t ) = S tot − (cid:88) i =3 x i ( t ) . ∂ z h i ( x ∗ , . . . , x ∗ ¯ d ) = − if i = 1 , − if i = 9 , and otherwise . Combining this with (31), we obtain the matrix − ∂ g ( x ∗ , . . . , x ∗ ¯ d ) · · · − ∂ ¯ d g ( x ∗ , . . . , x ∗ ¯ d ) − ∂ z h ( x ∗ , . . . , x ∗ ¯ d ) λ − ∂ g ( x ∗ , . . . , x ∗ ¯ d ) · · · − ∂ ¯ d g ( x ∗ , . . . , x ∗ ¯ d ) − ∂ z h ( x ∗ , . . . , x ∗ ¯ d )) ... − ∂ h ¯ d ( x ∗ , . . . , x ∗ ¯ d ) · · · λ − ∂ ¯ d g ¯ d ( x ∗ , . . . , x ∗ ¯ d ) − ∂ z h ¯ d ( x ∗ , . . . , x ∗ ¯ d )) shown in (26). By plugging in λ = 0 into the matrix and computing its determinant, we have H (0) = 3 . × . Since ¯ d = 9 , the third condition hold in Lemma 6.2.Consequently we show that all the conditions in Lemma 6.2 and show that Theorem hold 6.3,hence the linear stability of ( S , C , R , K ) follows. In this section, we use the expanded ACR controller (29) to control the 2-dimensional system A (cid:37) A + B −→ ∅ (cid:38) B, (32)which is introduced in Section 2.4 of the main text. The mass-action deterministic system associ-ated with this network is ddt a ( t ) = f ( a, b ) = − a ( t ) b ( t ) + 3 ,ddt b ( t ) = f ( a, b ) = − a ( t ) b ( t ) + 5 . Note that this system does not admit a positive steady state, as there does not exists ( a ∗ , b ∗ ) ∈ R > such that − a ∗ b ∗ = 0 and − a ∗ b ∗ = 0 . In particular, lim t →∞ ( b ( t ) − a ( t )) = ∞ since ddt ( b ( t ) − a ( t )) =2 . Furthermore the union system of (32) and the basic ACR system Z + A −→ Z and Z −→ A also does not admit a positive steady state. The associated mass-action system is ddt a ( t ) = − a ( t ) b ( t ) + 3 − z ( t )( a ( t ) − ,ddt b ( t ) = − a ( t ) b ( t ) + 5 ,ddt z ( t ) = z ( t )( a ( t ) − . ( a ∗ , b ∗ , z ∗ ) . By the last equation, a ∗ = 5 . Plugging a ∗ = 5 into a ( t ) in the first equation, we have b ∗ = . However, at this state, b is not stabilized as − a ∗ b ∗ + 5 = 2 , hence it contradicts to the assumption that ( a ∗ , b ∗ , z ∗ ) is a positive steady state.Now we consider the union system of (32) and the expanded ACR controller (29) introducedin Section 7. Z + A θ −→ Z, Z µ −→ A, Z + B α − (cid:42)(cid:41) − α Z, with general positive parameters κ , κ , κ , θ, µ, α and α . We use Lemma 7.1 to show this unionsystem admits a positive steady state under a mile condition. The associated mass-action system is ddt a ( t ) = − a ( t ) b ( t ) + 3 − z ( t )( θa ( t ) − µ ) ,ddt b ( t ) = − a ( t ) b ( t ) + 5 − z ( t )( α b ( t ) − α ) ,ddt z ( t ) = z ( t )( θa ( t ) − µ ) . It can be easily shown that for any positive constant c , there exists d = 3 d such that f ( c, d ) = 0 .Hence for a ∗ = µθ with arbitrary µ > and θ > , there exists b ∗ = 3 θµ such that f ( a ∗ , b ∗ ) = 0 .Since f ( a ∗ , b ∗ ) = 2 > , we tune the parameters α > and α > in (29) as they satisfy α b ∗ > α .Hence by Lemma 7.1, ( a ∗ , b ∗ , z ∗ ) = (cid:18) µθ , κ θκ µ , κ µ ( κ − κ ) α κ θ − α κ µ (cid:19) is a positive steady state. To control a stochastic system via an ACR controller, we rely on an approximation under multi-scaling model reduction as described in Section 2.4 and 2.5 of the main text. In this section weintroduce the formal procedures for generating a reduced model, and we introduce related theo-rems.
To formally define a reduced model of a given stochastic system, a notion of network projectionneeds to be introduced. The reduced system shown in Figure 4 b and the hybrid system shown in6 d are obtained through network projection. For example, consider the reaction network A κ − (cid:42)(cid:41) − κ B. (33)34uppose that for some parameters, species B is rarely produced or removed until time t = 1 . Inthis case, we approximate the distribution of A with the stochastic system associated with A κ B (0) −−−− (cid:42)(cid:41) −−−− κ . (34)Note that this reduced network is obtained by freezing the copy number of B at B (0) . We callnetwork (34) the projection of (33) by freezing species B at B (0) . As this example shows, networkprojection can be used to describe an asymptotic behavior of a subset of species.We define a projection function for complexes and reactions in ( S , C , R ) with S = S L ∪S H where S L = { S , S , . . . , S d } and S H = { S d +1 , S d +2 , . . . , S d + r } . In the later section, S L and S H would represent collections of species with low and high copy numbers, respectively.Let q L : Z d + r → Z d and q H : Z d + r → Z r be projection function such that for each v =( v , . . . , v d , v d +1 , . . . , v d + r ) T ∈ Z d + r , q L ( v ) = ( v , v , . . . , v d ) T ∈ Z d and q H ( v ) = ( v d +1 , v d +2 , . . . , v d + r ) T ∈ Z r . We use the projection function q L and q H for complexes and reaction. For example, for anetwork A + B → B with complexes A and B , we let S L = { A } . Then the complex A , B and thereaction A + B → B are represented with two dimensional vectors (1 , T , (0 , T and ( − , ,respectively. Then the projection of A , B and the reaction A + B → B are q L ((1 , T ) = 1 , q L ((0 , T ) = 0 and q L (( − , T ) = − which are associated with complexes A , and reaction A → , respectively. Hence by abusing notation, q L ( A ) = A, q L ( B ) = 0 and q L ( A + B → B ) = A → . Generally, we denote q L ( y ) the complex obtained by projection of the complexvector associated with a complex y . In this way, the projected network ( S L , C L , R L ) of the originalreaction network ( S , C , R ) by q L is defined to be S L = { S , . . . , S d } , C L = { q L ( y ) : y ∈ C} , and R L = { q L ( y ) → q L ( y (cid:48) ) : y → y (cid:48) ∈ R such that q L ( y (cid:48) ) − q L ( y ) (cid:54) = −→ } . (35)The rate constants of the projected network are defined with a given rate constants K of a givennetwork ( S , C , R ) . We inherit K to the projected network by incorporating the terms coming fromfreezing species S H at their initial count. For example, the rate constant of reaction A → in (34)is κ B (0) , where κ is inherited from reaction A + B → B .Hence for the set of reaction intensities K in a given network ( S , C , R , K ) , the set of reactionintensities for a projected network is K L = ¯ λ u ( x ) = (cid:88) y k → y (cid:48) k ∈R q L ( y k )=¯ y u ,q L ( y (cid:48) k )=¯ y (cid:48) u q H ( X (0)) q H ( y ) λ k ( x ) : ¯ y u → ¯ y (cid:48) u ∈ R L , (36)where u v = (cid:81) u v i i for the same dimensional non-negative vectors u and v and we use here the con-vention = 1 . The summation in (36) arises when multiple reactions in ( S , C , R ) are projectedinto the same reaction in R L . 35 .2 Stationary distribution approximation under multiscaling model reduc-tion The main idea of the stationary distribution approximation shown in Section 2.4 in the main textis the multiscaling model reduction. In this section, we introduce the multiscaling for a stochasticreaction network system with reaction propensities of constant order.Throughout this section we use the following notations. Let ( S , C , R , K ) be a network systemwith S = { S , S , . . . , S d + r } and let N be a scaling parameter. Let X N = ( X N , X N , . . . , X Nd + r ) be the associated scaled stochastic process with transition probabilities (12) such that each X Ni represents the counts of species S i . For a given initial condition X N (0) , let S L = { S i ∈ S : X Ni (0) = c i } and S H = { S i ∈ S : X Ni (0) = c i N } , (37)where c i ’s are positive constants.Then for a given collection of rate constants K , we scale the system with the following reactionintensities. K N = (cid:110) λ Nk = κ y → y (cid:48) N (cid:107) q H ( y ) (cid:107) λ k : λ k ∈ K (cid:111) , (38)where (cid:107) · (cid:107) is the 1-norm.For example, let A + B −→ be a reaction network system with S L = { A } and S H = { B } . Let x = ( x A , x B ) be a state of the associated stochastic process X N with a scaling parameter N . Thenthe for a given reaction intensity λ ( x ) = κx A x B of the reaction A + B → , we define a scaledreaction intensity λ N ( x ) = κN (cid:107) qH ( A + B ) (cid:107) λ ( x ) = κN x A x B . Then the scaled stochastic process X N has the transition probability P ( X N ( t + ∆ t ) = x + ( − , − T | X N ( t ) = x ) = λ N ( x )∆ t + h (∆) , where h is defined as (12).Having presented the above example, we now describe in more detail the procedure for a formalmultiscaling model reduction. Note that in the example above, as long as X B ( t ) is of order N , thepropensity is of order 1. Under this circumstance, we intuitively expect that X B would not besubstantially change because of the relatively low reaction propensity. Using this background, wecan approximate the distribution of a stochastic system through multiscaling model reduction. Theproof of the following theorem is provided in the separate paper [40]. Theorem 9.1.
Let X N = ( X N , X N , . . . , X Nd + r ) be the stochastic processes associated with ( S , C , R , K N ) where K N is as (38) . For an initial condition X (0) , suppose S = S L ∪ S H with S L and S H as in (37) . Let X be the associated stochastic process for the projected network system ( S L , C L , R L , K NL ) .Let further that p NL ( · , t ) and p ( · , t ) be the probability distributions for q L ( X N ) and X at time t ,respectively. Then for any A ⊂ Z d ≥ , we have | p N ( A, t ) − p ( A, t ) | = O ( N ν ) for some ν ∈ (0 , . (39) Remark 9.2.
In Theorem 9.1, if X admits a stationary distribution π , then | p N ( A, t ) − π ( A ) | ≤ | p N ( A, t ) − p ( A, t ) | + | p ( A, t ) − π ( A ) | . (40)Hence, for fixed t , if | p ( A, t ) − π ( A ) | is sufficiently small, then p NL ( A, t ) ≈ π ( A ) for large N .36or a special case, p NL ( A, t ) could be explicitly estimated with a Poisson distribution π Thefollowing corollary follows by Theorem 4.2 and Theorem 9.1.
Corollary 9.3.
Under the same conditions in Theorem 9.1, suppose that ( S L , C L , R L ) has zerodeficiency and is weakly reversible reaction network. Then for a positive steady state c ∈ R |S L | > ofthe deterministic counter part (11) of ( S L , C L , R L , K L ) , we have lim t →∞ lim N →∞ p NL ( x, t ) = M |S L | (cid:89) i =1 c x i i x i ! { x ∈ S } (41) where S is the state space, and M is the normalizing constant. In this section we consider the mean of the target species in a reduced system. As described in The-orem (4.2), when the reduced network ( S L , C L , R L ) has zero deficiency and is weakly reversible,the long-term behavior of the associated stochastic system follows a product form of Poissons (orconstrained Poissons). Moreover the mean of the system is determined by the positive steady stateof the deterministic counterpart as described in Theorem (4.2). Hence, we control the mean ofthe target species approximately by using the positive steady state value of the reduced system. Infact, for a certain reduced network obtained from an ACR system, the steady state value of an ACRspecies is preserved in the reduced network.Let π be a product form of Poisson distribution (or constrained Poissons) such as π ( x ) = M (cid:81) di =1 c xii x i ! { x ∈ S } defined on a state space S with some normalizing constant M and some c ∈ R d ≥ .Let π be the marginal distribution of the x -coordinate. If the support of π is whole Z ≥ then π is the Poisson distribution with rate c because π ( x ) = (cid:88) x ∈ Z ≥ · · · (cid:88) x d ∈ Z ≥ π ( x ) = 1 (cid:80) ∞ x i =0 c x x ! c x x ! = e − c c x x ! Thus the mean of the marginal distribution π is equal to c . Note that the reactions in the basicACR system (15) are always projected to µ − (cid:42)(cid:41) − θ X , the state space of X is always equal to Z ≥ in the projected network system. Hence the mean of the target species in the controlled system isclose to c , which is the positive steady state of the deterministic counter part.Similarly to the stability analysis carried out in Section 6, we assume a given network system ( S , C , R , K ) is ACR and admits conservation relations. With the same notation we used in Section6, let u i · x (0) = u i · x ( t ) = M i for all t (42)for some positive constants M i ’s and for some vectors u i ’s in R d> . Hence for each positive steadystate x ∗ of the system, if the i -th entry of x ∗ corresponds to a non-ACR species, then x ∗ i dependson the total concentrations M i . In this case, we denote x ∗ ( M ) a positive steady state on the stoi-chiometry class S x (0) , where M = ( M , M , . . . , M k ) such that u i · x (0) = M i . We also denotethe entry of x ∗ ( M ) by x ∗ i = x ∗ i ( M ) . With these notations, the following theorem states that if areduced network is obtained by freezing some non-ACR species of an ACR system, then the steadystate value of an ACR species is preserved in the reduced system.37 heorem 9.4. Let ( S , C , R , K ) be a mass-action ACR system with x (cid:48) ( t ) = f ( x ( t )) . Suppose ( S , C , R , K ) admits the conservation laws (42) . We split S = { S , . . . , S d , S d +1 , . . . , S d + r } intotwo disjoint subsets S L = { S , . . . , S d } and S H = { S d +1 , . . . , S d + r } , where none of the species in S H is an ACR species. Let { x ∗ ( M ) : M = ( M , M , . . . , M k ) ∈ R k> } be the family of positivesteady states of the system ( S , C , R , K ) . For the dynamics x ( t ) of ( S , C , R , K ) , suppose there exista (cid:102) M = ( (cid:102) M , . . . , (cid:102) M k ) such that x d + i (0) = x ∗ d + i ( (cid:102) M , . . . , (cid:102) M k ) for i = 1 , , . . . , r. (43) Then the projected system ( S L , C L , R L , K L ) with the initial condition q L ( x (0)) has a positivesteady state value at ¯ x ∗ = q L ( x ∗ ( (cid:102) M )) . In particular, if S i ∈ S L is an ACR species with theACR value x ∗ i , then ¯ x ∗ i = x ∗ i .Proof. Let ¯ x (cid:48) ( t ) = ¯ f (¯ x ( t )) be the deterministic system of ( S L , C L , R L , K L ) . Note that ( S L , C L , R L ) is obtained by freezing the species in S H from ( S , C , R ) . Recall by the definition of K L (38) thateach element K L is a summation of some rate constants in K multiplied by the initial values ofspecies in S H . Hence for any positive x i ’s, we have ¯ f i ( x , . . . , x d ) = f i ( x , . . . , x d , x d +1 (0) , . . . , x d + r (0)) , i = 1 , , . . . , d. For the (cid:102) M , we have f ( x ∗ ) = f ( x ∗ ( (cid:102) M ) , . . . , x ∗ d , x ∗ d +1 ( (cid:102) M ) , . . . , x ∗ d + r ( (cid:102) M ))= f ( x ∗ ( (cid:102) M ) , . . . , x ∗ d ( (cid:102) M ) , x d +1 (0) , . . . , x d + r (0))= ¯ f ( x ∗ ( (cid:102) M ) , . . . , x ∗ d ( (cid:102) M )) . Hence ¯ x ∗ = ( x ∗ ( (cid:102) M ) , . . . , x ∗ d ( (cid:102) M )) = q L ( x ∗ ( (cid:102) M )) is a positive steady state of ( S L , C L , R L , K L ) .Suppose S i ∈ S L is an ACR species with the ACR value x ∗ i . Since the ACR value is in-dependent of any conservative quantities M , M , . . . , M k by definition, we have x ∗ i ( (cid:102) M ) = x ∗ i .Consequently, ¯ x ∗ i = x ∗ i . Example 9.1.
Consider a network ( S , C , R , K ) C + A D + AA B Z + A Z Z A θ µ
There are two conservation relations a (0)+ b (0)+ z (0) = a ( t )+ b ( t )+ z ( t ) = M and c (0)+ d (0) = c ( t ) + d ( t ) = M . The mass action deterministic dynamics associated with ( S , C , R , K ) is a (cid:48) ( t ) = − a ( t ) + b ( t ) − z ( t )( θa ( t ) − µ ) ,b (cid:48) ( t ) = a ( t ) − b ( t ) ,z (cid:48) ( t ) = z ( t )( θa ( t ) − µ ) ,c (cid:48) ( t ) = − c ( t ) a ( t ) + d ( t ) a ( t ) ,d (cid:48) ( t ) = − d ( t ) a ( t ) + c ( t ) a ( t ) , ence the positive steady state x ∗ = ( a ∗ , b ∗ , c ∗ , d ∗ , z ∗ ) = ( µθ , µθ , M , M , M − µθ ) . This impliesthat species A and B are ACR species.We split S into S L = { A, B, C } and S H = { D, Z } . For a given initial condition x (0) =( a (0) , b (0) , c (0) , d (0) , z (0)) , we have the following reduced network ( S L , C L , R L , K L )0 A BC + Aµz (0) θz (0) 11 d (0) 1 With M = 2 µθ + z (0) and M = 2 d (0) , the condition (43) holds for species E and Z . Henceby Theorem 9.4, the positive steady state values of A and B in ( S L , C L , R L , K L ) are same as thepositive steady state values in ( S , C , R , K ) , which are µθ for both species. (cid:52) In this section, we show that the basic ACR controller (15) can also be used to control a stochasticsystem under the classical scaling regime. To show that the target species in the scaled stochasticsystem follows approximately a Poisson distribution, we use a hybrid type system. This frameworkis introduced in [20].We define γ k = 1 if q H ( y k ) > for a reaction y k → y (cid:48) k ∈ R , otherwise γ k = 0 . As shown inSection 9.2, for a given set of intensity functions K , we define new scaled reaction intensities as ¯ K N = (cid:110) ¯ λ Nk = κ y → y (cid:48) N (cid:107) q H ( y ) (cid:107) − γ k λ k : λ k ∈ K (cid:111) . (44)Suppose S = S L ∪ S H as (37) for a reaction network ( S , C , R ) . Suppose further that the pro-jected network ( S L , C L , R L ) has deficiency of zero and is weakly reversible. According to [20,Theorem 3.8 and Corollary 4.4], the system dynamics for ( S , C , R , ˜ K N ) can be approximatelystudied with a hybrid system. The copy numbers of the species in S L are modeled with a stochas-tic process, and the distribution of q L ( X N ) is approximated with a product form of Poissons ata finite time t . On the other hand, the concentration of the species in S H follows the determin-istic dynamics. These two processes are coupled as the rate constants of the stochastic part aredetermined by the deterministic part, and the kinetics of the deterministic part also depends on thestationary mean of the stochastic part. We denote by ( S H , C H , R H , ¯ K H ) the deterministic part ofthe hybrid system. More precise definition of ( S H , C H , R H , ¯ K H ) can be found in [20].In the scaled stochastic system modeled with ¯ K N defined at (44), let m N ( t ) be the mean copynumber of the target species X at time t such that lim N →∞ m N ( t ) = m ( t ) . That is, m ( t ) approximatesthe mean of X in the scaled stochastic system. The following Lemma provides the conditionsguaranteeing that m ( t ) converges in time to the desired value in a controlled system with the basicACR controller. As we assumed for the basic ACR system in Section 5, we suppose that X is ina given reaction network and Z is a newly introduced control species. Lemma 9.5.
For a given reaction network, let ( S , C , R ) be the union of the given network and thebasic ACR system (15) . Let X N be the stochastic process associated with ( S , C , R , ¯ K N ) , where K N is a set of scaled reaction intensities defined as (44) . For an initial state of X N (0) , suppose X ∈ S L and Z ∈ S H . Suppose further that the conditions of Theorem 3.8 in [20] hold. If theconcentration of Z in the deterministic part ( S H , C H , R H , ¯ K H ) converges to a positive steady state,as t → ∞ , then lim t →∞ m ( t ) = µθ .Proof. Since Z is newly introduced species with the basic ACR system, Z is regulated with onlytwo reactions ∅ ← Z → A in ( S H , C H , R H , ¯ K H ) . Let z ( t ) be the concentration of species Z attime t . According to Theorem 3.8 in [20], the rate of reaction ∅ ← Z is µz ( t ) , and the rate ofreaction Z → A is θm ( t ) z ( t ) . Then z ( t ) solves a differential equation, ddt z ( t ) = z ( t ) ( θm ( t ) − µ ) , Hence if lim t →∞ z ( t ) = z ∗ for some z ∗ ∈ (0 , ∞ ) , then lim t →∞ ( θm ( t ) − µ ) = 0 .In Section 10.2, we use this lemma to show that the mean of the species X in the hybrid systemin Figure (6)d converges to µθ . By equation (39), it is important to show the term | p ( A, t ) − π ( A ) | is small for the stationary dis-tribution approximation (41) with a Poisson distribution. Basically the term | p ( A, t ) − π ( A ) | tendsto zero as t → ∞ by the ergodic theorem [51]. In this section we introduce one of the most well-known theoretical frameworks, the so-called Foster-Lyapunov criterion [52] for the convergenceof | p N ( A, t ) − π ( A ) | in t . The following theorem is a version of Theorem 6.1 in [52]. Theorem 9.6 (Foster-Lyapunov criterion for exponential ergodicity) . Let X be a continuous-timeMarkov chain on a countable state space S with the infinitesimal generator A (14) . Suppose thereexists a positive function V on S satisfying the following conditions.1. V ( x ) → ∞ as | x | → ∞ , and2. There are positive constants a and b such that A V ( x ) ≤ − aV ( x ) + b for all x ∈ S . (45) Then X admits a unique stationary distribution π . Furthermore, there exists positive constants η and C such that for any measurable set A and any state x , | P ( X ( t ) ∈ A | X (0) = x ) − π ( A ) | ≤ C ( V ( x ) + 1) e − ηt . For a finite time t , if the projected network system in Theorem 9.1 satisfies the conditions inTheorem 9.6, then the term | p ( A, t ) − π ( A ) | in the right hand side of (40) can be small. Therefore | p N ( A, t ) − π ( A ) | in (40) can eventually be small enough with sufficiently large N .40 In this section we provide the initial reaction propensities in the receptor-ligand signaling modeldescribed in Figure 4a. Note that we model the reaction propensities with stochastic mass-actionkinetics (13). Letting N = 10 , we set the initial values L (0) = 1 . N, P (0) = 100 , R (0) =2 , Z (0) = N , and we set zero initial values for the rest of species. The parameters are κ N =1 . / (1 . N ) , κ N = 1 . , κ N = 1 . , κ N = 1 . , κ N = 1 . , κ N = 1 . , κ N = 1 . , κ N =1 . , κ N = 1 . , κ N = 1 . , κ N = 1 . , κ N = 0 . , θ = 1 /N and µ = 5 /N . Then for the initialcondition X (0) = ( L (0) , R (0) , R (0) , D (0) , D (0) , D (0) , D (0) , Z (0)) , we have the followingpropensities of reaction L + R → R, R → L + R and R → D . λ L + R → R ( X (0)) = κ L (0) R (0) = 2 . , λ R → L + R ( X (0)) = κ R (0) = 2 . λ R → D ( X (0)) = κ R (0)( R (0) −
1) = 2 . . The propensities of all remaining reactions are zero since the initial counts of
D, D , D and D are zero. As we mentioned in Section 2.4 of the main text, the initial propensities of each reactionis relatively small to the initial counts of species L and Z . Hence it needs a longer time to substan-tially deviate the counts of L and Z . Over a short term interval, each of L and Z in the associatedstochastic process behaves as a constant function (Figure 4c). Hence the reduced model obtainedby freezing L and Z at their initial values approximates the original controlled system.By using Theorem 9.1, we show this approximation more precisely. We let ( S , C , R ) be thecontrolled system in Figure 4 a. We choose the initial values and the set K N of parameters asintroduced above with the scaling parameter N . Note that K N satisfies (38). We split S into S L = { R , R, D, D , D , D } and S H = { L, Z } . Then the reduced network ( S L , C L , R L , K L ) isthe network in Figure 4 b. The deficiency of the reduced network is as the number of complexes is , there are two linkage classes, and the rank of the stoichiometry matrix is . The reduced networkis also weakly reversible because each linkage class is strongly connected. Hence Corollary 9.3implies the distribution p NL of S L species at t = 150 in ( S , C , R , K ) is estimated by a product formof Poissons as described at (41).In order to approximate the mean of R in the controlled system ( S , C , R ) , we show that thepositive steady state value of R in ( S L , C L , R L , K L ) is µθ = 5 . Theorem 9.4 can be used to showthat the mean, but using deficiency zero condition of the reduced network provides a much simplerway. Since ( S L , C L , R L , K L ) has zero deficiency and is weakly reversible, the associated massaction deterministic dynamics is complex balanced [16, 17]. This means that for each complex, all‘in-flows‘ and ‘out-flows‘ are balanced at each positive steady state. Hence for the zero complexin ( S L , C L , R L , K L ) , the in-flow is θr ∗ and the out-flow is µ for the positive steady value r ∗ of R .Therefore they are balanced at r ∗ = µθ .We now investigate the accuracy of this approximation. Let p NL ( · , t ) and p ( · , t ) denote the dis-tribution of species S L in the controlled receptor-ligand system and the distribution of the reducedsystem, respectively. We further let π be the product form of Poissons stationary distribution ofthe reduced system. As shown in (40), for a small error between p NL ( · , t ) and π with a fixed time t = 150 , we need a fast convergence for p ( · , t ) to π .With a slight modification on the reduced network, we show that how to use the Foster-Lyapunov criterion in Theorem 9.6 to show the convergence rate of p ( · , t ) to π in time t . In order41o construct a Lyapunov function explicitly, we add a degradation of D to the reduced model inFigure 4b. Hence let ( S L , C L , R L , K L ) be a system described with the following reaction network. R R R DD D D κ κ κ κ κ κ κ κ κ κ µθκ (46)Let x = ( x , x , x , x , x , x ) be a vector each of whose entries represents the copy numbersof R , R, D, D , D and D , respectively. We use a linear Lyapunov function V ( x ) = (cid:80) i =1 v i x i with some positive vector v = ( v , v , . . . , v d ) . The work in [53] exploited details about linearLyapunov functions for stochastic reaction networks. By the definition of the generator A (14), wehave A V ( x ) = κ ( v − v ) x + ( κ ( v − v ) + κ ( v − v )) x + ( κ ( v − v ) + κ ( v − v )) x + ( κ ( v − v ) + κ (2 v − v ) − κ v ) x + + κ ( v − v ) x + ( κ ( v − v ) − θv ) x + ( κ ( v − v )) x + µ. (47)Let h i = v i +1 − v i for i = 1 , , , and h = v − v . We will choose h i ’s with which all thecoefficients on the right hand side of (47) are strictly negative numbers. Hence, we choose h i ’ssuch that h < , h > , κ κ h − κ κ h < v = 2 v + 2 h + h κ κ h < h , κ κ h < h , and h > . (48)With a sufficiently large v , we can find h i ’s satisfying (48). Then we have A V ( x ) = − c x − c (cid:48) x + c x − c x − c x − c x − c x , for some c i > for i = 1 , , . . . , and c (cid:48) > . Hence (45) holds, and for the distribution p ( · , t ) of (46) and its stationary distribution π , we have the exponential decay of | p ( A, t ) − π ( A ) | for any A ⊂ Z ≥ , as t → ∞ . Recall that we used the hybrid system shown in Figure 6d to approximate the distribution of X inthe controlled dimer-catalyzer system that we introduced in Figure 6a. In this section, by usingLemma 9.5 we prove that m ( t ) , the mean of X at time t in the hybrid system, converges to µθ undera mild condition, as t → ∞ . 42et x ( t ) , c ( t ) , c p ( t ) , c pp ( t ) , x ∗ ( t ) and z ( t ) be the solution to the deterministic part of the hybridsystem shown in Figure 6d. Let further m ( t ) = κ x ∗ ( t ) + µz ( t ) κ c pp ( t ) + θz ( t ) and (cid:96) ( t ) = κ c ( t ) κ x ( t ) be the meanof X and X at time t , respectively. Then the deterministic system is governed by the followingsystem of ordinary differential equations, ddt x ( t ) = κ c ( t ) − κ x ( t ) (cid:96) ( t ) = 0 ,ddt c ( t ) = κ x ( t ) (cid:96) ( t ) + κ c p ( t ) − ( κ + κ ) c ( t ) = κ c p ( t ) − κ c ( t ) ,ddt c p ( t ) = κ c ( t ) + κ c pp ( t ) − ( κ + κ ) c p ( t ) ,ddt c pp ( t ) = κ c p ( t ) − κ c pp ( t ) ,ddt x ∗ ( t ) = κ c pp ( t ) m ( t ) − κ x ∗ ( t ) , and ddt z ( t ) = z ( t )( θm ( t ) − µ ) . (49)Let x (0) , c (0) , c p (0) , c pp (0) , x ∗ (0) and z (0) be the initial values of the system (49). Note that ddt x ∗ ( t ) + ddt z ( t ) = 0 and ddt c ( t ) + ddt c p ( t ) + ddt c pp ( t ) = 0 for all t . We let M = x ∗ (0) + z (0) = x ( t ) + z ( t ) and L = c (0) + c p ( t ) = c ( t ) + c p ( t ) + c pp ( t ) be the conserved quantities.By using the second, third and fourth equations in (49), there is a single positive steady state, (¯ c, ¯ c p , ¯ c pp ) = (cid:18) κ κ Lκ κ + κ κ + κ κ , κ κ κ Lκ κ + κ κ κ + κ κ κ , κ κ κ κ Lκ κ + κ κ κ + κ κ κ κ (cid:19) , for x ( t ) , c ( t ) , c p ( t ) and c pp ( t ) . In order to study the stability of this positive steady state, weconsider the following linear system for c, c p and c pp , ddt cc p c pp = − κ κ κ − κ − κ κ κ − κ cc p c pp . The eigenvalues of the matrix above is λ = 0 ,λ = −
12 ( κ + κ + κ + κ ) − (cid:112) ( κ + κ − κ − κ ) + 4 κ κ ,λ = −
12 ( κ + κ + κ + κ ) + 12 (cid:112) ( κ + κ − κ − κ ) + 4 κ κ . Since ( κ + κ + κ + κ ) > ( κ + κ − κ − κ ) + 4 κ κ , the eigenvalues λ and λ are strictlynegative. Thus c ( t ) c p ( t ) c pp ( t ) = v + v e − λ t + v e − λ t , v i ’s are the corresponding eigenvectors. In particular, v = (¯ c, ¯ c p , ¯ c pp ) . This implies that lim t →∞ c pp ( t ) = ¯ c pp .Now in the following proposition we introduce a sufficient condition for z ( t ) to converge tosome positive steady state. Proposition 10.1.
Suppose that θκ M > µκ ¯ c pp . Then z ( t ) in (49) converges to θκ M − µκ ¯ c pp θκ ,as t → ∞ .Proof. By using the conservative quantity M = x ∗ ( t ) + z ( t ) , the differential equation for z in (49)can be written as ddt z ( t ) = z ( t ) (cid:18) θκ M − µκ c pp ( t ) − θκ z ( t ) κ c pp ( t ) + θz ( t ) (cid:19) = θκ z ( t ) κ c pp ( t ) + θz ( t ) ( α ( t ) − z ( t )) , where α ( t ) = M − µκ θκ c pp ( t ) .Let ¯ α = M − µκ θκ ¯ c pp . Let also (cid:15) > be an arbitrarily small number.Since lim t →∞ c pp ( t ) = ¯ c pp ,there exists a T > such that | α ( t ) − ¯ α | < (cid:15) when t > T . We denote R − (cid:15) = ¯ α − (cid:15) and R + (cid:15) = ¯ α + 2 (cid:15) . Then clearly α ( t ) ∈ [ R − (cid:15) , R + (cid:15) ] for all t > T .In the rest of the proof, we suppose t > T . Note that if z ( t ) < R − (cid:15) , then ddt z ( t ) = θκ z ( t ) κ c pp ( t ) + θz ( t ) ( α ( t ) − z ( t )) > θκ z ( t ) κ c pp ( t ) + θz ( t ) ( α ( t ) − ¯ α + 2 (cid:15) ) > θκ z ( t ) κ L + θM (cid:15) Note further that if z ( t ) > R + (cid:15) , then ddt z ( t ) = θκ z ( t ) κ c pp ( t ) + θz ( t ) ( α ( t ) − z ( t )) < θκ z ( t ) κ c pp ( t ) + θz ( t ) ( α ( t ) − ¯ α − (cid:15) ) < − θκ z ( t ) κ L + θM (cid:15) Therefore there exists T = inf { t ≥ T : z ( t ) ∈ [ R − (cid:15) , R + (cid:15) ] } < ∞ . Let T = inf { t > T : z ( t ) = ¯ α } . We first consider the case z ( T ) = R − (cid:15) and T < ∞ . If there exists t > T such that z ( t ) > α ( t ) , then by using the continuity of z ( t ) , z ( t ) = z ( T ) + (cid:90) tT dds z ( s ) ds = α ( t ) + (cid:90) tT θκ z ( s ) κ c pp ( s ) + θz ( s ) ( α ( s ) − z ( s )) ds < α ( t ) . Hence it is contradiction to z ( t ) > α ( t ) . Thus R − epsilon < z ( t ) < α ( t ) for all t > T . Now weconsider the case z ( T ) = R − (cid:15) and T = ∞ . Then for any t > T , it follows that z ( t ) = z ( T ) + (cid:90) tT dds z ( s ) ds = R − epsilon + (cid:90) tT θκ z ( s ) κ c pp ( s ) + θz ( s ) ( α ( s ) − z ( s )) ds > R − (cid:15) as z ( t ) < α ( t ) for any t > T . Hence R − (cid:15) < z ( t ) < α ( t ) for all t > T . Hence we concludethat when z ( T ) = R − (cid:15) , there exists T (cid:48) ≥ T ≥ T such that such that α ( t ) < z ( t ) < R + (cid:15) for all t > T (cid:48) . In the same way, it follows that when z ( T ) = R + (cid:15) there exists T (cid:48) ≥ T ≥ T such that α ( t ) < z ( t ) < R + (cid:15) for all t > T .Consequently, there exists T (cid:48) ≥ T such that z ( t ) ∈ [ R − (cid:15) , R + (cid:15) ] for all t > T (cid:48) . Since (cid:15) can bechosen arbitrarily small with sufficiently large T , z ( t ) converges to ¯ α , as t → ∞ . Obviously ¯ α ispositive by the assumption θκ M > µκ ¯ c pp . 44roposition 10.1 shows that for any choice of system parameters κ i , if either the initial value of x ∗ or z is large enough, then z ( t ) converges to a positive steady state, as t goes to ∞ . Hence by us-ing Lemma 9.5, we conclude that lim t →∞ m ( t ) = µθ if we input sufficiently large initial concentrationof the control species Z . References [1] David F. Anderson, Gheorghe Craciun, and Thomas G. Kurtz. Product-form stationary distri-butions for deficiency zero chemical reaction networks.
Bull. Math. Biol. , 72(8):1947–1970,2010.[2] German Enciso and Jinsu Kim. Embracing noise in chemical reaction networks.
Bulletin ofmathematical biology , 81(5):1261–1267, 2019.[3] Hye-Won Kang, Wasiur R. KhudaBukhsh, Heinz Koeppl, and Grzegorz A. Rempala. Quasi-steady-state approximations derived from the stochastic model of enzyme kinetics.
Bull.Math. Bio. , 2019.[4] Johan Paulsson. Summing up the noise in gene networks.
Nature , 427:415–418, 2004.[5] Guy Shinar and Martin Feinberg. Structural sources of robustness in biochemical reactionnetworks.
Science , 327(5971):1389–1391, 2010.[6] Martin Feinberg. Lectures on Chemical Reaction networks. https://crnt.osu.edu/LecturesOnReactionNetworks , 1979.[7] Corentin Briat, Ankit Gupta, and Mustafa Khammash. Antithetic integral feedback ensuresrobust perfect adaptation in noisy biomolecular networks.
Cell systems , 2(1):15–26, 2016.[8] Stephanie K Aoki, Gabriele Lillacci, Ankit Gupta, Armin Baumschlager, David Schwein-gruber, and Mustafa Khammash. A universal biomolecular integral feedback controller forrobust perfect adaptation.
Nature , page 1, 2019.[9] Corentin Briat, Ankit Gupta, and Mustafa Khammash. Antithetic proportional-integral feed-back for reduced variance and improved control performance of stochastic reaction networks.
J. R. Soc. Interface , 15:20180079, 2018.[10] Adam Arkin, John Ross, and Harley H. McAdams. Stochastic kinetic analysis of develop-mental pathway bifurcation in phage lambda-infected Escherichia coli cells.
Genetics , 149:1633–1648, 1998.[11] A Becskei, Benjamin B. Kaufmann, and A Van Oudenaarden. Contributions of low moleculenumber and chromosomal positioning to stochastic gene expression.
Nature Genetics , 37(9):937–944, 2005.[12] Michael B. Elowitz, Arnold J. Levin, Eric D. Siggia, and Peter S. Swain. Stochastic GeneExpression in a Single Cell.
Science , 297(5584):1183–1186, 2002.4513] D Huh and Johan Paulsson. Non-genetic heterogeneity from stochastic partitioning at celldivision.
J. Nat. Genet. , 43(2):95–100, 2011.[14] H´edia Maamar, Arjun Raj, and David Dubnau. Noise in Gene Expression Determines CellFate in Bacillus subtilis.
Science , 317(5837):526–529, 2007.[15] S Uphoff, N D Lord, L Potvin-Trottier, B Okumus, D J Sherratt, and J Paulsson. Stochasticactivation of a DNA damage response causes cell-to-cell mutation rate variation.
Science ,351(6277):1094–1097, 2016.[16] Friedrich Josef Maria Horn. Necessary and sufficient conditions for complex balancing inchemical kinetics.
Arch. Rat. Mech. Anal. , 49(3):172–186, 1972.[17] M Feinberg. Complex balancing in general kinetic systems.
Arch. Rational Mech. Anal. , 49:187–194, 1972.[18] Gheorghe Craciun. Toric differential inclusions and a proof of the global attractor conjecture. arXiv preprint arXiv:1501.02860 , 2015.[19] German A. Enciso. Transient absolute robustness in stochastic biochemical networks.
J. R.Soc. Interface , 13:20160475, 2016.[20] David F Anderson, Daniele Cappelletti, and Thomas G Kurtz. Finite time distributions ofstochastically modeled chemical systems with absolute concentration robustness.
SIAM Jour-nal on Applied Dynamical Systems , 16(3):1309–1339, 2017.[21] David F Anderson, Germ´an A Enciso, and Matthew D Johnston. Stochastic analysis of bio-chemical reaction networks with absolute concentration robustness.
Royal Society Interface ,11:20130943, 2014.[22] Thomas G Kurtz. The Relationship between Stochastic and Deterministic Models for Chem-ical Reactions.
J. Chem. Phys. , 57(7):2976–2978, 1972.[23] Karen Ball, Thomas G Kurtz, Lea Popovic, and Greg Rempala. Asymptotic analysis ofmultiscale approximations to reaction networks.
The Annals of Applied Probability , 16(4):1925–1961, 2006.[24] Hye-Won Kang and Thomas G Kurtz. Separation of time-scales and model reduction forstochastic reaction networks.
Annals of Applied Probability , 23(2):529–583, 2013.[25] Boris Y Rubinstein, Henry H Mattingly, Alexander M Berezhkovskii, and Stanislav YShvartsman. Long-term dynamics of multisite phosphorylation.
Molecular biology of thecell , 27(14):2331–2340, 2016.[26] Carsten Conradi and Anne Shiu. A global convergence result for processive multisite phos-phorylation systems.
Bulletin of mathematical biology , 77(1):126–155, 2015.[27] Nida Obatake, Anne Shiu, Xiaoxian Tang, and Angelica Torres. Oscillations and bistabilityin a model of erk regulation. arXiv preprint arXiv:1903.02617 , 2019.4628] Eric Batchelor and Mark Goulian. Robustness and the cycle of phosphorylation and dephos-phorylation in a two-component regulatory system.
Proceedings of the National Academy ofSciences , 100(2):691–696, 2003.[29] Nicola Normanno, Antonella De Luca, Caterina Bianco, Luigi Strizzi, Mario Mancino, Mon-ica R Maiello, Adele Carotenuto, Gianfranco De Feo, Francesco Caponigro, and David SSalomon. Epidermal growth factor receptor (EGFR) signaling in cancer.
Gene , 366(1):2–16,2006.[30] Sarah R Needham, Selene K Roberts, Anton Arkhipov, Venkatesh P Mysore, Christopher JTynan, Laura C Zanetti-Domingues, Eric T Kim, Valeria Losasso, Dimitrios Korovesis,Michael Hirsch, et al. EGFR oligomerization organizes kinase-active dimers into compe-tent signalling platforms.
Nature communications , 7:13307, 2016.[31] James E Ferrell Jr. Perfect and near-perfect adaptation in cell signaling.
Cell Systems , 2(2):62–67, 2016.[32] Wenzhe Ma, Ala Trusina, Hana El-Samad, Wendell A Lim, and Chao Tang. Defining networktopologies that can achieve biochemical adaptation.
Cell , 138(4):760–773, 2009.[33] Fangzhou Xiao and John C Doyle. Robust perfect adaptation in biomolecular reaction net-works. In , pages 4345–4352. IEEE,2018. doi: 10.1109/CDC.2018.8619101.[34] Tau-Mu Yi, Yun Huang, Melvin I Simon, and John Doyle. Robust perfect adaptation inbacterial chemotaxis through integral feedback control.
Proceedings of the National Academyof Sciences , 97(9):4649–4653, 2000.[35] Ankit Gupta and Mustafa Khammash. A finite state projection algorithm for the stationarysolution of the chemical master equation.
The Journal of Chemical Physics , 147(15), 2017.[36] Jae Kyoung Kim, Grzegorz A Rempala, and Hye-Won Kang. Reduction for stochastic bio-chemical reaction networks with multiscale conservations.
Multiscale Modeling & Simula-tion , 15(4):1376–1403, 2017.[37] Jae Kyoung Kim and Eduardo D Sontag. Reduction of multiscale stochastic biochemi-cal reaction networks using exact moment derivation.
PLoS computational biology , 13(6):e1005571, 2017.[38] Brian Munsky and Khammash Mustafa. The finite state projection approach for the analysisof stochastic noise in gene networks Dissertation.
IEEETrans. Autom. Contr , 53:201–214,2008.[39] Dan T. Gillespie. Exact stochastic simulation of coupled chemical reactions.
J. Phys. Chem. ,81(25):2340–2361, 1977.[40] German Enciso and Jinsu Kim. Constant Order Multiscaling Reduction for Stochastic Reac-tion Networks. arXiv , 2019. 4741] Luigi Preziosi. Hybrid and multiscale modelling.
Journal of mathematical biology , 53(6):977–978, 2006.[42] Hao Ge, Hong Qian, and X Sunney Xie. Stochastic phenotype transition of a single cell in anintermediate region of gene state switching.
Physical review letters , 114(7):078101, 2015.[43] Yen Ting Lin and Nicolas E Buchler. Efficient analysis of stochastic gene dynamics in thenon-adiabatic regime using piecewise deterministic markov processes.
Journal of The RoyalSociety Interface , 15(138):20170804, 2018.[44] Hao Ge, Pingping Wu, Hong Qian, and Sunney Xiaoliang Xie. Relatively slow stochasticgene-state switching in the presence of positive feedback significantly broadens the region ofbimodality through stabilizing the uninduced phenotypic state.
PLoS computational biology ,14(3):e1006051, 2018.[45] Colin S Gillespie. Moment-closure approximations for mass-action models.
IET systemsbiology , 3(1):52–58, 2009.[46] Daniele Cappelletti, Ankit Gupta, and Mustafa Khammash. A hidden integral structure en-dows absolute concentration robust systems with resilience to dynamical concentration dis-turbances. arXiv preprint arXiv:1910.05531 , 2019.[47] Oren Shoval, Lea Goentoro, Yuval Hart, Avi Mayo, Eduardo Sontag, and Uri Alon. Fold-change detection and scalar symmetry of sensory input fields.
Proceedings of the NationalAcademy of Sciences , 107(36):15995–16000, 2010.[48] Guy Shinar, Ron Milo, Mar´ıa Rodr´ıguez Mart´ınez, and Uri Alon. Input–output robustness insimple bacterial signaling systems.
Proceedings of the National Academy of Sciences , 104(50):19931–19935, 2007.[49] Joseph P Dexter and Jeremy Gunawardena. Dimerization and bifunctionality confer robust-ness to the isocitrate dehydrogenase regulatory system in escherichia coli.
Journal of Biolog-ical Chemistry , 288(8):5770–5778, 2013.[50] Stewart N Ethier and Thomas G Kurtz.
Markov Processes: Characterization and Conver-gence . John Wiley & Sons, New York, 1986.[51] James Norris.
Markov Chains . Cambridge University Press, 1997.[52] Sean P. Meyn and Richard L. Tweedie. Stability of Markovian Processes III : Foster-Lyapunov Criteria for Continuous-Time Processes.
Advances in Applied Probability , 25(3):518–548, 1993. ISSN 00018678. doi: 10.2307/1427522. URL .[53] Ankit Gupta, Corentin Briat, and Mustafa Khammash. A scalable computational frameworkfor establishing long-term behavior of stochastic reaction networks.