Modular Composition of Gene Transcription Networks

Predicting the dynamic behavior of a large network from that of the composing modules is a central problem in systems and synthetic biology. Yet, this predictive ability is still largely missing because modules display context-dependent behavior. One cause of context-dependence is retroactivity, a phenomenon similar to loading that influences in non-trivial ways the dynamic performance of a module upon connection to other modules. Here, we establish an analysis framework for gene transcription networks that explicitly accounts for retroactivity. Specifically, a module’s key properties are encoded by three retroactivity matrices: internal, scaling, and mixing retroactivity. All of them have a physical interpretation and can be computed from macroscopic parameters (dissociation constants and promoter concentrations) and from the modules’ topology. The internal retroactivity quantifies the effect of intramodular connections on an isolated module’s dynamics. The scaling and mixing retroactivity establish how intermodular connections change the dynamics of connected modules. Based on these matrices and on the dynamics of modules in isolation, we can accurately predict how loading will affect the behavior of an arbitrary interconnection of modules. We illustrate implications of internal, scaling, and mixing retroactivity on the performance of recurrent network motifs, including negative autoregulation, combinatorial regulation, two-gene clocks, the toggle switch, and the single-input motif. We further provide a quantitative metric that determines how robust the dynamic behavior of a module is to interconnection with other modules. This metric can be employed both to evaluate the extent of modularity of natural networks and to establish concrete design guidelines to minimize retroactivity between modules in synthetic systems. Citation: Gyorgy A, Del Vecchio D (2014) Modular Composition of Gene Transcription Networks. PLoS Comput Biol 10(3): e1003486. doi:10.1371/ journal.pcbi.1003486 Editor: Stanislav Shvartsman, Princeton University, United States of America Received May 17, 2013; Accepted January 10, 2014; Published March 13, 2014 Copyright: 2014 Gyorgy, Del Vecchio. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Funding: Grant #FA9550-12-1-0129. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The authors have declared that no competing interests exist. * E-mail: ddv@mit.edu


Introduction
The ability to accurately predict the behavior of a complex system from that of the composing modules has been instrumental to the development of engineering systems.It has been proposed that biological networks may have a modular organization similar to that of engineered systems and that core processes, or motifs, have been conserved through the course of evolution and across different contexts [1], [2], [3], [4], [5].In addition to having profound consequences from an evolutionary perspective, this view implies that biology can be understood, just like engineering, in a modular fashion [6].To predict the behavior of a network from that of its composing modules, it is certainly desirable that the salient properties of modules do not change upon connection with other modules.This modularity property is especially important in a bottom-up approach to engineer biological systems, in which small systems are combined to create larger ones [7], [8].
Unfortunately, despite the fact that biological networks are rich of frequently repeated motifs, suggesting a modular organization, a module's behavior is often affected by its context [9].Contextdependence is due to a number of different factors.These include unknown regulatory interactions between the module and its surrounding systems; various effects that the module has on the cell network, such as metabolic burden [10], effects on cell growth [11], and competition for shared resources [12]; and loading effects associated with known regulatory linkages between the module and the surrounding systems, a phenomenon known as retroactivity [13], [14].As a result, our current ability of predicting the emergent behavior of a network from that of the composing modules remains limited.This inability is a central problem in systems biology and especially daunting for synthetic biology, in which circuits need to be re-designed through a lengthy and ad hoc process every time they are inserted in a different context [15].
In the phenomenon known as retroactivity, a downstream module perturbs the dynamic state of its upstream module in the process of receiving information from the latter [13], [14].These effects are due to the fact that, upon interconnection, a species of the upstream module becomes temporarily unavailable for the reactions that make up the upstream module, changing the upstream module's dynamics.The resulting perturbations can have dramatic effects on the upstream module's behavior.For example, in experiments in gene circuits in Escherichia coli, a few fold ratio in gene copy number between the upstream module and the downstream target results in more than 40% change in the upstream module's response time [16].More intriguing effects take place when the upstream module is a complex dynamical system such as an oscillator.In particular, experiments in transcriptional circuits in vitro showed that the frequency and amplitude of a clock's oscillations can be largely affected by a load [17] and computational studies on the genetic activator-repressor clock of [18] further revealed that just a few additional targets for the activator impose enough load to quench oscillations.Surprisingly, adding a few targets for the repressor can restore the stable limit cycle [19].Retroactivity has also been experimentally demonstrated in signaling networks in vitro [20] and in the MAPK cascade in vivo [21].In particular, it was shown in [19] that a few fold ratio between the amounts of the upstream and downstream system's proteins can lead to more than triple the response time of the upstream system.
In this paper, we provide a quantitative framework to accurately predict how and the extent to which retroactivity will change a module's temporal dynamics for general gene transcription networks and illustrate the implications on a number of recurrent network motifs.We demonstrate that the dynamic effects of loading due to interconnections can be fully captured by three retroactivity matrices.The first is the internal retroactivity, which accounts for loading due to intramodular connections.We illustrate that due to internal retroactivity, negative autoregulation can surprisingly slow down the temporal response of a gene as opposed to speeding it up, as previously reported [22]; perturbations applied at one node can lead to a response at another node even in the absence of a regulatory path from the first node to the second, having consequences relevant for network identification techniques (e.g., reviewed in [23]); and an oscillator design can fail even in the presence of small retroactivity.The other two matrices, which we call scaling and mixing retroactivity, account for loading due to intermodular connections.We illustrate that because of the scaling retroactivity, the switching characteristics of a genetic toggle switch can be substantially affected when the toggle switch is inserted in a multi-module system such as that proposed for artificial tissue homeostasis in [24].The interplay between scaling and internal retroactivity plays a role in performance/robustness trade-offs, which we illustrate considering the single-input motif [5].Using these retroactivities, we further provide a metric establishing the robustness of a module's behavior to interconnection.This metric can be explicitly calculated as a function of measurable biochemical parameters, and it can be used both for evaluating the extent of modularity of natural networks and for designing synthetic circuits modularly.
Our work is complementary to but different from studies focusing on partitioning large transcription networks into modules using graph-theoretic approaches [13], [25], [26].Instead, our main objective is to develop a general framework to accurately predict both the quantitative and the qualitative behavior of interconnected modules from their behavior in isolation and from key physical properties (internal, scaling, and mixing retroactivity).In this sense, our approach is closer to that of disciplines in biochemical systems analysis, such as metabolic control analysis (MCA) [27], [28].However, while MCA is primarily focused on steady state and near-equilibrium behavior, our approach considers global nonlinear dynamics evolving possibly far from equilibrium situations.This paper is organized as follows.We first introduce a general mechanistic model for gene transcription networks to explain the physical origin of retroactivity and to formulate the main question of the paper (System Model and Problem Formulation).We then provide the two main results of the paper (Results).These are obtained by reducing the mechanistic model through the use of time scale separation (leading to models of the same dimension as those based on Hill functions), in which only macroscopic parameters and protein concentrations appear.In these reduced models, the retroactivity matrices naturally arise, whose practical implications are illustrated on five different application examples.

System Model and Problem Formulation
We begin by introducing a standard mechanistic model for gene transcription networks, which includes protein production, decay, and reversible binding reactions between transcription factors (TFs) and promoter sites, required for transcriptional regulation.Specifically, transcription networks are usually viewed as the input/output interconnection of fundamental building blocks called transcriptional components.A transcriptional component takes a number of TFs as inputs, and produces a single TF as an output.The input TFs form complexes with promoter sites in the transcriptional component through reversible binding reactions to regulate the production of the output TF, through the process of gene expression (for details, see Methods).To simplify the notation, we treat gene expression as a one-step process, neglecting mRNA dynamics.This assumption is based on the fact that mRNA dynamics occur on a time scale much faster than protein production/decay [1].In addition to this, including mRNA dynamics is not relevant for the study of retroactivity, and would yield only minor changes in our results (see Methods).
Within a transcription network, we identify a transcriptional component with a node.Consequently, a transcription network is a set of interconnected nodes in which node x i represents the transcriptional component producing TF x i .There is a directed edge from node x j to x i if x j is a TF regulating the activity of the promoter controlling the expression of x i [29], in which case we call x j a parent of x i .Activation and repression are denoted by ? and a, respectively.Modules are a set of connected nodes.Modules communicate with each other by having TFs produced in one module regulate the expression of TFs produced in a different module.When a node x i is inside the module, we call the corresponding TF x i an internal TF, while when node x i is outside the module we call the corresponding TF x i an external TF.Further, we identify external TFs that are parents to internal TFs as inputs to the module.Let x, u and c denote the concentration vector of internal TFs, inputs and TF-promoter complexes, respectively.According to [30], we can write the dynamics of the module as

Author Summary
Biological modules are inherently context-dependent as the input/output behavior of a module often changes upon connection with other modules.One source of context-dependence is retroactivity, a loading phenomenon by which a downstream system affects the behavior of an upstream system upon interconnection.This fact renders it difficult to predict how modules will behave once connected to each other.In this paper, we propose a general modeling framework for gene transcription networks to accurately predict how retroactivity affects the dynamic behavior of interconnected modules, based on salient physical properties of the same modules in isolation.We illustrate how our framework predicts surprising and counter-intuitive dynamic properties of naturally occurring network structures, which cannot be captured by existing models of the same dimension.We describe implications of our findings on the bottom-up approach to designing synthetic circuits, and on the topdown approach to identifying functional modules in natural networks, revealing trade-offs between robustness to interconnection and dynamic performance.Our framework carries substantial conceptual analogies with electrical network theory based on equivalent representations.We believe that the framework we have proposed, also based on equivalent network representations, can be similarly useful for the analysis and design of biological networks.
where N st is the stoichiometry matrix and v is the reaction flux vector.The reactions are either protein production/decay or binding/unbinding reactions.Therefore, we partition v into r Ã and r, representing the reaction flux vectors corresponding to production/decay and binding/unbinding reactions, respectively (see Methods).We assume that the DNA copy number is conserved, therefore, we can rewrite (1) as where the upper left block matrix in N st is the zero matrix as DNA is not produced/degraded.As a result, with g x,c which we call the isolated dynamics of a module.Next, consider the case when the module is inserted into a network, which we call the context of the module.We represent all the quantities related to the context with an overbar.Let x x and c c denote the concentration vector of TFs and promoter complexes of the context, respectively.Furthermore, denote by r r Ã and r r the reaction flux vectors corresponding to production/decay and binding/unbinding reactions between TFs and promoters in the context of the module, respectively.Then, the dynamics of the species in the module (c and x) and in the context ( c c and x x) can be written as where the upper left block matrix is zero as DNA is assumed to be a conserved species.Furthermore, since r and r r encapsulate the binding/unbinding reactions in the module and in its context, respectively, the off-diagonal block matrices in the upper right block matrix are zero.Similarly, as r Ã and r r Ã encapsulate the production/decay reactions in the module and its context, respectively, the off-diagonal block matrices in the lower left block matrix are zero.Finally, the stoichiometry matrix E E represents how internal TFs of the module participate in binding/unbinding reactions in the context of the module (E can be interpreted similarly).
With s~ E E r r describing the effective rate of change of x due to intermodular binding reactions, we obtain which we call the connected dynamics of a module.We refer to s as the retroactivity to the output of the module, encompassing retroactivity applied to the module due to the context of the module.Similarly, we call r the retroactivity to the input of a module, representing retroactivity originating inside the module.The general interconnection of a group of modules can be treated similarly (Figure 1).As an example of the implications of retroactivity s on the module's dynamic behavior, consider Figure 2.For the purpose of illustration, assume that f 1 t ð Þ and f f 1 t ð Þ, external inputs to x 1 and x x 1 (see Methods), are periodic (in general, they can be arbitrary Figure 1.The dynamics of a module depend on the module's context.Downstream modules change the dynamics of an upstream module by applying a load.The effect of this load is captured by the retroactivity to the output s of the upstream module, which is the weighted sum of the retroactivity to the input r (i) of the downstream modules.doi:10.1371/journal.pcbi.1003486.g001time-varying signals).When the module is not connected to its context (Figure 2A), its output is periodic (Figure 2B).Upon interconnection with its context (Figure 2C), due to the retroactivity to the output s applied by the context, the output of the module changes significantly (Figure 2D).Hence, connection with the context leads to a dramatic departure of the dynamics of the module from its behavior in isolation.This example illustrates that retroactivity s significantly alters the dynamic behavior of modules after interconnection, therefore, it cannot be neglected if accurate prediction of temporal dynamics is required.Unfortunately, model (4) provides little analytical insight into how measurable parameters and interconnection topology affect retroactivity.
The aim of this paper is to provide a model that captures the effects of retroactivity, unlike standard regulatory network models of the same dimension based on Hill functions [1].Specifically, we seek a model that explicitly describes the change in the dynamics of a module once it is arbitrarily connected to other modules in the network.This model is only a function of measurable biochemical parameters, TF concentrations, and interconnection topology.

Results
We first characterize the effect of intramodular connections on an isolated module's dynamics.We then analytically quantify the effects of intermodular connections on a module's behavior.Finally, we determine a metric of robustness to interconnection quantifying the extent by which the dynamics of a module are affected by its context.We demonstrate the use of our framework and its implications on network motifs taken from the literature.
The main technical assumptions in what follows are that (a) there is a separation of time scale between production/degradation of proteins and the reversible binding reactions between TFs and DNA, and that (b) the corresponding quasi-steady state is locally exponentially stable.Assumption (a) is justified by the fact that gene expression is on the time scale of minutes to hours while binding reactions are on the second to subsecond time scale [3].Assumption (b) is implicitly made any time Hill function-based models are used in gene regulatory networks.In addition to these technical assumptions, to simplify notation, we model gene expression as a one-step process, however, a more detailed description of transcription/translation would not yield any changes to the main results (see Methods).

Effect of Intramodular Connections
Here, we focus on a single module without inputs and describe how retroactivity among nodes, modeled by Br in (2), affects the module's dynamics.To this end, we provide a model that well approximates the isolated module dynamics, in which only measurable macroscopic parameters appear, such as dissociation constants and TF concentrations.We then present implications of this model for negative autoregulation, combinatorial regulation and the activator-repressor clock of [18].
Employing assumptions (a)-(b), we obtain the first main result of the paper as follows.Let x~(x 1 ,x 2 , . . .,x N ) T denote the vector of concentrations of internal TFs, then the dynamics well approximate the dynamics of x in (2) in the isolated module with h(x)~f where f i represents external perturbations to x i (inducer, noise, or disturbance, f i (t):0 unless specified otherwise), d i is the decay rate of x i , and H i (p i ) is the Hill function modeling the production rate of x i , regulated by the parents p i of x i .We call R i (p i ) the retroactivity of node x i .For the most common binding types, Figure 3 shows the expressions of H i (p i ) and R i (p i ) (for their definition, see Methods).The binary matrix V i has as many columns as the number of nodes in the module, and as many rows as the number of parents of x i , such that its (j,k) element is 1 if the j th parent of x i is x k , otherwise the entry is zero.That is, an entry in the following matrix x 1 x 2 . . .
. . .If node x i has no parents, its node retroactivity is not defined.In the single parent case, node x i has one parent, y binding as an n-multimer with dissociation constant k y .In the case of independent, competitive and cooperative binding, node x i has two parents, y and z, binding as multimers with multimerization factors n and m, respectively, together with dissociation constants k y and k z , respectively.The total concentration of the promoter of x i is denoted by g i .The production rates p i,0 , p i,1 , p i,2 and p i,3 correspond to the promoter complexes without parents, with y only, with z only, and with both y and z, respectively.For details, see Supporting Text S3. doi:10.1371/journal.pcbi.1003486.g003 is 1 if the species indexing the corresponding row and column are the same, otherwise the entry is zero, yielding p i ~Vi x.Finally, W is the set of nodes having parents from inside the module.For the derivation of this result, see Theorem 1 in Supporting Text S2.
We call R the internal retroactivity of the module as it describes how retroactivity among the nodes internal to the module affects the isolated module dynamics.When R~0, we have _ x x~h(x), the commonly used Hill function-based model for gene transcription networks [3].It is possible to show that h(x) represents the rate of change of total (free and bound) TFs (see Supporting Text S2).Hence, (6) describes how changes in the total concentration of TFs h(x) relate to changes _ x x in the concentration of free TFs.Specifically, to change the concentration of free TFs by one unit, the module has to change the total concentration of TFs by (IzR) units, as R units are ''spent on'' changing the concentration of bound TFs.Having R~0 implies that the module's effort on affecting the total concentration of TFs is entirely spent on changing the concentration of free TFs.By contrast, R k k?? implies that no matter how much the total concentration of TFs changes, it is not possible to achieve any changes in the free concentration of some of the TFs.Therefore, the internal retroactivity R describes how ''stiff'' the module is against changes in x due to loading applied by internal connections.
The entries of R i (p i ) have the following physical interpretation.Consider first a module with the autoregulated node x 1 , that is, x 1 has a single parent: itself.The retroactivity of node x 1 is R 1 (x 1 )~a, where a is given in Figure 3.In this case, we obtain R(x 1 )~a by ( 6), so that (5) yields _ x x 1 ~1 1za h(x 1 ).Hence, the greater a, the harder to change the concentration of free x 1 by changing its total concentration (the ''stiffer'' the node), and the temporal dynamics of x 1 become slower.The retroactivity R 1 of node x 1 can be increased by increasing its DNA copy number g 1 or by decreasing the dissociation constant k 1 of x 1 .For a node with two parents, we provide the explicit formula for R i in Figure 3 in the case of the most frequent binding types, so that here we simply write The diagonal entries b and e in ( 7) can be interpreted similarly to a, while the off-diagonal entries can be interpreted as follows.Having cw0 means that the second parent facilitates the binding of the first, whereas cv0 represents blockage (d can be interpreted similarly with the parents having reverse roles).Therefore, we have c~d~0 in the case of independent binding (Figure 3), as the parents bind to different sites.By contrast, we have c,dƒ0 in the case of competitive binding (Figure 3), since the parents are competing for the same binding sites, forcing each other to unbind.Following a similar reasoning, we obtain c,d §0 in the case of cooperative binding (Figure 3).Notice that R i is scaled by the total concentration of promoter g i , which can be changed, for example, in synthetic circuits by changing the plasmid copy number.

Practical Implications of Intramodular Connections
In order to illustrate the effects of intramodular connections, we consider three recurrent network motifs in gene transcription networks: (i) negative autoregulation of a gene, (ii) combinatorial regulation of a gene by two TFs, and (iii) the activator-repressor clock of [18].
Negative autoregulation.One of the most frequent network motifs in gene transcription networks is negative autoregulation, as over 40% of known Escherichia coli TFs are autorepressed [29].Earlier studies concluded that negative autoregulation makes the response of a gene faster [22].Here, we demonstrate that in the case of significant retroactivity, negative autoregulation can actually slow down the response of a gene.To this end, consider a module consisting of the single node x 1 , and analyze first the case when its production is constitutive with promoter concentration g 1 , production rate constant p 1,0 and decay rate d 1 .Then, the dynamics of x 1 are given by _ x x 1 ~g1 p 1,0 {d 1 x 1 .In the case of negative autoregulation, x 1 has itself as the only parent.Let k 1 denote the dissociation constant of x 1 and assume it binds as a monomer repressing its own production (so that n~1 and p 1,1 ~0 in Figure 3).According to Figure 3, we have and W~fx 1 g, yielding from ( 6) R(x 1 )~R 1 (x 1 ) and h(x 1 )~H 1 (x 1 ) {d 1 x 1 , so that (5) results in This expression indicates that negative autoregulation yields two changes in the dynamics.First, protein production changes from g 1 p 1,0 to the Hill function H 1 (x 1 ).Second, the dynamics are The red and blue plots denote the cases with and without negative autoregulation, respectively, whereas the green plot represents the case of negative autoregulation neglecting retroactivity (R~0 in (8)).Simulation parameters are d 1 ~1hr {1 , k 1 ~10nM, together with p 1,0 ~20hr {1 , p 1,0 ~10hr {1 , p 1,0 ~1hr {1 for A, B and C, respectively.
To carry out a meaningful comparison between the unregulated and regulated systems, we compare the response time of systems with the same steady state.To do so, we pick the same value of g 1 in the case of the regulated systems (g 1 ~15nM, g 1 ~30nM, g 1 ~300nM for A, B and C, respectively), but a different one for for the unregulated system (g 1 ~2:5nM, g 1 ~5nM, g 1 ~50nM for A, B and C, respectively), such that the steady states match (see Methods for parameter ranges).Decreasing p 1,0 (lower production rate constant) while increasing g 1 (higher DNA copy number) results in slower response, as internal retroactivity increases.doi:10.1371/journal.pcbi.1003486.g004 premultiplied by 1zR ð Þ {1 , which is the effect of internal retroactivity.
As it was shown in [22], the response time of the regulated system without retroactivity is smaller than that of the unregulated system.When considering internal retroactivity, however, the response time increases, as the absolute value of _ x x 1 decreases with increased R according to (8).Specifically, the response time with R is greater than without R since Rw0.That is, while the Hill function makes the response faster, internal retroactivity has an antagonistic effect, so that negative autoregulation can render the response slower than that of the unregulated system, as illustrated in Figure 4. Furthermore, if p 1,0 g 1 is kept constant, the response time of both the unregulated (blue) and the regulated system without retroactivity (green) remain the same, together with the steady states.By contrast, increasing g 1 (and decreasing p 1,0 ) makes the internal retroactivity R greater (since R is proportional to g 1 ), while the contribution of the Hill function remains unchanged.As a result, the response of the regulated system with retroactivity (red) becomes slower as we increase g 1 (and decrease p 1,0 ).This is illustrated in Figure 4 with different (g 1 ,p 1,0 ) pairs.Note that p 1,0 can be decreased, for example, by decreasing the ribosome binding site (RBS) strength, whereas g 1 can be increased by increasing the gene copy number.
Combinatorial regulation.As a second example, we consider a single gene co-regulated by two TFs (Figure 5A).This topology appears in recurrent network motifs, such as the feedforward-loop, the bi-fan and the dense overlapping regulon [5].Here, we show that a perturbation introduced in one of the parents (blue in Figure 5A) can affect the concentration of the other parent (red node in Figure 5A), even in the absence of a regulatory path between the two.
Referring to ( 5)-( 6), note that x 3 is the only node with parents (W~fx 3 g), so that R(x)~V T  3 R 3 V 3 .Using (6) with where the entries of R 3 are given in Figure 3 (depending on the binding type at x 3 ) together with H 3 (x 1 ,x 2 ), the dynamics in ( 5) take the form : This expression implies that unless d~0, a perturbation f 1 (Figure 5B) in x 1 yields a subsequent perturbation in x 2 .In the case of independent binding, we have d~0, and as a result, no perturbation is observed in x 2 .In the case of competitive binding, instead, we have dv0, so that perturbations f 1 in x 1 yield perturbations of the same sign in x 2 , that is, x 1 acts as if it were an activator of x 2 (Figure 5C).In the case of cooperative binding, instead, we have dw0.As a result, perturbations in x 1 yield perturbations in x 2 of opposite sign (Figure 5D), which implies that x 1 behaves as if it were a repressor of x 2 .As d is proportional to g 3 (Figure 3), higher DNA copy number for x 3 yields greater pulses in x 2 subsequent to an equal perturbation in x 1 .Interestingly, if we view x 2 as the output of the module, the module has the adaptation property with respect to its input x 1 (or f 1 ).That is, retroactivity enables to respond to sudden changes in input stimuli, while adapting to constant stimulus values.
Activator-repressor clock.One common clock design is based on two TFs, one of which is an activator and the other is a repressor [18], [31], [32].Here, we illustrate the effect of internal retroactivity on the functioning of the clock design of [18] depicted in Figure 6A.In particular, x 1 activates the production of both TFs, whereas x 2 represses the production of x 1 through competitive binding.Consequently, the network topology is captured by the binary matrices V 1 ~I and V 2 ~1 0 ½ , whereas h(x) and R(x) can be constructed by considering (H 1 (x 1 ,x 2 ), H 2 (x 1 )) and (R 1 (x 1 ,x 2 ), R 2 (x 1 )), respectively, in Figure 3. Here, we write R 2 ~a, while the entries of R 1 are denoted by b, c, d and e, as in (7).Then, we obtain that (5) takes the form (C) In the case of competitive binding, increasing the concentration of free x 1 yields more of x 1 bound to the promoter of x 3 , forcing some of the molecules of x 2 to unbind, thus increasing the free concentration x 2 .Consequently, x 1 acts as if it were an activator of x 2 .(D) By contrast, in the case of cooperative binding, when the binding of x 1 must precede that of x 2 , pulses in x 1 yield pulses of the opposite sign in x 2 .Consequently, x 1 acts as if it were a repressor of x 2 .Simulation parameters are: g 1 ~g2 ~10nM, g 3 ~20nM, d 1 ~d2 ~1hr {1 , p 1,0 ~0, p 2,0 ~10hr {1 , k 1 ~k2 ~1nM 4 , and both x 1 and x 2 bind as tetramers.doi:10.1371/journal.pcbi.1003486.g005 It was previously shown [33] that the principle for the clock to oscillate is that the activator dynamics are sufficiently faster than the repressor dynamics (so that the unique equilibrium point is unstable).Equation (9) shows that the activator and repressor dynamics are slowed down asymmetrically (diagonal terms in ½IzR(x) {1 ), and that they become coupled (off-diagonal terms in ½IzR(x) {1 , c,d=0), because of internal retroactivity.In particular, in the case when c,d%1ze%1zazb, the activator would slow down compared to the repressor.Based on the principle of functioning of the clock, we should expect that this could stabilize the equilibrium point, quenching the oscillations as a consequence.In fact, oscillations disappear even if the circuit is assembled on DNA with a single copy (g 1 ~g2 ~2nM), as it can be observed in Figure 6D.Therefore, accounting for internal retroactivity is particularly important in synthetic biology during the design process when circuit parameters and parts are chosen for obtaining the desired behavior.An effective way to restore the limit cycle in the clock, yielding sustained oscillations, is to render the repressor dynamics slower with respect to the activator dynamics.This can be obtained by adding extra DNA binding sites for the repressor [19], as shown in Figure 6B.In fact, in this case, we have R 3 (x 2 )w0 given in Figure 3, which, due to (5), will yield the following change in (9): instead of e, we will have ezR 3 we, rendering the dynamics of the repressor slower with respect to the activator dynamics.As a result, the equilibrium point becomes unstable, restoring the limit cycle, verified by simulation in Figure 6E.Further studies on specific systems have investigated the effects of TF/promoter binding on the dynamics of loop oscillators, such as the repressilator [34].

Effect of Intermodular Connections
After investigating how retroactivity due to intramodular connections affect a single module's dynamics, we next determine how the dynamics of a module change when the module is inserted into its context.To this end, we first extend the model in (5) to the case in which the module has external TFs as inputs.Hence, let u~(u 1 ,u 2 , . . .,u W ) T denote the concentration vector of TFs external to the module.With this, we obtain that the dynamics well approximate the dynamics of x in (2) with where V is the set of nodes having parents from outside the module (external TFs), and the binary matrix D i has as many columns as the number of inputs of the module, and as many rows as the number of parents of x i , such that its (j,k) element is 1 if the j th parent of x i is u k , otherwise the entry is zero.That is, an entry in the following matrix u 1 u 2 . . .
. . . is 1 if the species indexing the corresponding row and column are the same, otherwise the entry is zero, yielding p i ~½ V i D i ( x T u T ) T .Furthermore, note that in the presence of input u, both h(:) and R(:) given in (6) depend on x and u, as some of the parents of internal TFs are external TFs.For the derivation of this result, see Theorem 2 in Supporting Text S2.
Before stating the main result of this section, we first provide the interpretation of Q. Recall that h(x,u)~0 implies that the total concentrations of internal TFs are constant.In this case, (10) , where x is the concentration vector of free internal TFs.This means that the concentrations of free internal TFs can still be changed subsequent to changes in the external TFs (input), despite the fact that the total concentration (free and bound) of internal TFs remains unaffected.Therefore, Q Figure 6.Neglecting internal retroactivity could falsely predict that the activator-repressor clock will display sustained oscillations.(A) The module consists of the activator protein x 1 (dimer) and the repressor protein x 2 (tetramer), with dissociation constants k 1 and k 2 , respectively.(B) An extra node x 3 is introduced as a target for the repressor.(C) Without accounting for internal retroactivity, the module in A exhibits sustained oscillations.(D) When internal retroactivity is included for the module in A, however, the equilibrium point is stabilized and the limit cycle disappears.(E) Oscillations can be restored by applying a load on the repressor (module in B) with concentration g, so that the repressor dynamics are slowed down.Simulation parameters: p 1,0 ~0:04hr {1 , p 1,1 ~2500hr {1 , p 1,2 ~0hr {1 , p 2,0 ~0:004hr {1 , p 2,1 ~100hr {1 , d 1 ~1hr {1 , d 2 ~0:6hr {1 , g 1 ~g2 ~2nM, k 1 ~1nM 2 , k 2 ~1nM 4 and g~40nM.doi:10.1371/journal.pcbi.1003486.g006captures the phenomenon by which external TFs force internal TFs to bind/unbind, for instance, by competing for the same binding sites.Having Q~0 means that external TFs do not affect the binding/unbinding of internal TFs, which is the case, for example, when all bindings are independent.Thus, we call Q the external retroactivity of the module.
The main result of this section describes how the context of a module affects the module's dynamics due to retroactivity.Specifically, we consider the module of interest and we represent the rest of the network, the module's context, as a different module.As previously, we use the overbar to denote that a quantity belongs to the context.With this, we obtain that the dynamics ! |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} isolated dynamics of the module and of its context ð12Þ well approximate the dynamics of x and x x in (4) in the module connected to the context with : where N and N N denote the number of nodes in the module and in the context, respectively.Furthermore, the binary matrix U has as many rows as the number of inputs of the module, and as many columns as the number of nodes in the context, such that its (j,k) element is 1 if the j th input of the module is the k th internal TF of the context (u j ~ x x k ), otherwise the entry is zero.That is, an entry in the following matrix x x 1 x x 2 . . . . . . is 1 if the species indexing the corresponding row and column are the same, otherwise the entry is zero, yielding u~U x x.The quantities corresponding to the context, that is, S S, M M and U U are defined similarly with the only difference that variables with and without overbar have to be swapped (for instance, N and N N have to be swapped in (13)).For the derivation of this result, see Theorem 3 in Supporting Text S2.
We next provide the interpretation of the scaling and mixing retroactivity.The reduced order model (12) describes how retroactivity between the module and the context affects each other's dynamics.Note that zero matrices S, M, S S and M M lead to no alteration in the dynamics upon interconnection.To further deepen the implications of these matrices and their physical meaning, note that when M M~0, the dynamics of the module after interconnection become that is, S S determines how the isolated dynamics of the module get ''scaled'' upon interconnection.Therefore, we call S S the scaling retroactivity of the context, accounting for the loading that the context applies on the module as some of the TFs of the module are taken up by promoter complexes in its context (we obtain S S~0 if nodes in the context do not have parents in the module, that is, if V V~1).Since the dynamics of the context enter into the module's dynamics through M M, we call M M the mixing retroactivity of the context, referring to the ''mixing'' of the dynamics of the module and that of its context.The mixing retroactivity M M establishes how internal TFs force external TFs to bind/unbind, so that M M~0 can be ensured if the binding of parents from the module is independent from that of the parents from the context.This holds if nodes in the context are not allowed to have parents in both the module and in the context ( V V\ W W~1). When M M=0, a perturbation applied in the context can result in a response in the upstream module, even without TFs in the context regulating TFs in the module, leading to a counter-intuitive transmission of signals from downstream (context) to upstream (module).
With this, we can explain the simulation results in Figure 2D by analyzing S S and M M for the system in Figure 2C.Let R 1 ~a and let R R 2 be defined as in (7), where a, b, c, d and e are given in Figure 3.Then, we have R~a by (13) and S S~b and M M~( c 0 ) by ( 6).Hence, expression (12) yields describing the dynamics of x 1 upon interconnection with its context, where _ x x 1 ~f (x) and _ x x x x 1 ~ f f 1 ( x x,x, _ x x) describe the dynamics of x 1 and x x 1 , respectively, when the module and the context are not connected to each other.If M M~0, then (15) reduces to that is, the context rescales the dynamics of the module.The smaller (1za)=(1zazb), that is, the greater the scaling retroactivity b~ S S, the greater the effect of this scaling.Note that since the scaling factor is smaller than 1 (unless the scaling retroactivity is zero, i.e., b~0), the effect of the scaling retroactivity of the context in this case is to make the temporal dynamics of the module slower.Once M M=0, in addition to this sclaing effect, the dynamics of the context appear in the dynamics of the module (Figure 2D).Referring to (15), we can quantify the effect of the context on the module, considering the ratio c=(1za).The greater c=(1za), that is, the greater M M, the stronger the contribution of the context compared to that of the module to the dynamics of the module upon interconnection.Here, both b and c increase, for instance, with the copy number of x x 2 (Figure 3).
Connecting the module to its context such that x 1 and x x 1 are competing for the same binding sites is less desirable than employing independent binding, as the dynamics of the context (downstream system) can suppress the dynamics of the module (upstream system).Dismantling the dynamics of the module will ''misinform'' other downstream systems in the network that are regulated by x 1 .From a design perspective, multi-module systems should be designed and analyzed such that the modules have zero mixing retroactivity.This can be achieved, for instance, by avoiding non-independent binding at the interface nodes (at x x 2 in Figure 2B).However, since completely independent binding can be hard to realize in the case of combinatorial regulation, nodes integrating signals from different modules should not be placed into the input layer (nodes having parents from other modules), yielding V V\ W W~1.This can be achieved by introducing an extra input node in the downstream module (see Supporting Figure S1).
Next, we quantify the difference between the isolated and connected module behavior.In particular, we provide a metric of the change in the dynamics of a module upon interconnection with its context, dependent on R and S S and under the assumption that M M~0 (obtained, for instance, by avoiding mixing parents from the module and the context).The isolated dynamics of the module can be well approximated by the reduced order model _ x x~f (x,u, _ u u) in (10), and let x(t) denote its solution.Once we connect the module to its context, its dynamics change according to ( 14), which we write as _ x x~f f (x,u, _ u u) and let x x(t) denote the corresponding solution.Using the sub-multiplicative property of the induced 2-norm, we have that the percentage change of the dynamics upon interconnection can be bounded from above as follows: Furthermore, with m m §0 independent of x and u, such that m(x,u)ƒm m, we obtain that that is, m m provides an upper bound on the percentage change in the dynamics of the module, and on the difference in the trajectories of the module upon interconnection with its context.Furthermore, by using the properties of the induced 2-norm, we obtain that we can pick provided that s max ( S S)vs min (IzR) for all x and x x, where s min (IzR) and s max ( S S) denote the smallest singular value of (IzR) and largest singular value of S S, respectively.For the mathematical derivations, see Supporting Text S2.This suggests that the module becomes more robust to interconnection by increasing min x, x x s min (IzR) or by decreasing max x, x x s max ( S S).Such a metric can be used both in the analysis and in the design of complex gene transcription networks as follows.Given any network and a desired module size N (number of nodes within the module), we can identify the module that has the least value of m m, that is, the module with the greatest guaranteed robustness to interconnection.Furthermore, we can also evaluate existing partitionings based on other measures (e.g., edge betweenness [35], its extension to directed graphs with nonuniform weights [36], round trip distance [25] or retroactivity [37]) with respect to robustness to interconnection.From a design point of view, one can design multi-module systems such that internal, scaling and mixing retroactivities allow for low values of m m, leading to modules that behave almost the same when connected or isolated.

Practical Implications of Intermodular Connections
We next illustrate the effect of intermodular connections on the dynamics of interconnected modules, considering both a synthetic genetic module that is being employed in a number of applications and a natural recurring network motif.
Toggle switch.Here, we consider the toggle switch of [38], a bistable system that can be permanently switched between two steady states upon presentation of a transient input perturbation.This module has been proposed for synthetic biology applications in biosensing (see, for example, [34], [39]).In this paper, we consider the toggle switch inserted into the context of the synthetic circuit for controlling tissue homeostasis as proposed in [24], and investigate how the context of the toggle affects its switching characteristics.Figure 7A illustrates the toggle switch in isolation, whereas Figure 7B shows the configuration when connected to the context [24].Note that all nodes, both in the toggle switch and in its context, have a single parent.Therefore, ) are given in Figure 3.We first consider the model of the toggle switch when not connected to its context (Figure 7A).Since the toggle switch has no input, its isolated dynamics are described by (5), where Next, we consider the toggle switch connected to its context (Figure 7B).As nodes in the toggle switch have no parents from outside it, we have S~0 and M~0 by (13).Nodes in the context have no parents in the context, leading to R R~0 from ( 6), and to Q Q~0, referring to (11).With this, the isolated dynamics of the context are given by according to (10).The fact that nodes in the context do not have mixed parents from the toggle switch and from the context results in M M~0 from (13).With ) §0, the dynamics of the toggle switch once connected to the context (Figure 7B) are given by according to (12), so that the dynamics of x 1 and x 2 are unaffected if a~0 and b~0, respectively.When a,bw0, both the x 1 and x 2 dynamics become slower upon interconnection, so that the response to an input stimulation will also be slower.As a consequence, upon removal of the stimulation, the displacement in the toggle state may not be sufficient to trigger a switch.This is illustrated in Figure 7C-D.In order to recover the switch, a wider pulse is required (Figure 7E) to compensate for the slow-down due to the context (also, note that the switching dynamics are slower than in the isolated case).As a result, even if the toggle had been characterized in isolation, it would fail to function as expected when inserted into its context.Note that we have m ), where a represents the amount of load on x 1 imposed by the context compared to that by the module, and b can be interpreted similarly.The greater a (or b), the slower the dynamics of x 1 (of x 2 ) become upon interconnection with the context.Greater a and b yield greater m m, suggesting decreased robustness to interconnection, verified by the simulation results.
Single-input motif.As a second example, we focus on a recurrent motif in gene transcription networks, called the singleinput motif [5].The single-input motif is defined by a set of operons (context) controlled by a single TF (module), which is usually autoregulated (Figure 8A).It is found in a number of instances, including the temporal program controlling protein assembly in the flagella biosynthesis [40].Here, we show that the dynamic performance (speed) of the module and its robustness to interconnection with its context are not independent, and that this trade-off can be analyzed by focusing on the interplay between the internal retroactivity R of the module and the scaling retroactivity S S of the context.
The isolated dynamics of the module are given in ( 8), which we write here as _ x x 1 ~f (x 1 ).Furthermore, we have D D i ~1 for i~1,2, . . .,l and U U~1, so that S S(x 1 )~P l i~1 R R i (x 1 ) by ( 13), where l is the number of nodes in the context and R R i (x 1 ) is given in Figure 3 (single parent).Consequently, upon interconnection, the dynamics of the module change according to (14) as where m(x 1 )~ S S(x 1 )=½1zR(x 1 )z S S(x 1 ) and equals the expression in (16).The smaller m(x 1 ), the more robust the module to interconnection.Note that R(x 1 ) is proportional to g 1 , therefore, while increasing R(x 1 ) makes the module slower (Figure 8B), it also makes it more robust to interconnection (Figure 8C).It was previously shown that negative autoregulation increases robustness to perturbations [41].Here, we have further shown that increasing the internal retroactivity R of the module provides an additional mechanism to increase robustness to interconnection, at the price of slower response.For a fixed steady state (the product of p 1,0 and g 1 is held constant), smaller p 1,0 yields greater g 1 , that is, increased R and, in turn, smaller m.From a design perspective, if speed is a priority, one should choose a strong RBS with a low copy number plasmid, or alternatively, a promoter with high dissociation constant k 1 .By contrast, if robustness to interconnection is central, a weak RBS with a high copy number plasmid (or with low k 1 ) is a better choice.If both speed and robustness to interconnection are desired, other design approaches may be required, such as the incorporation of insulator devices, as proposed in other works [42].The toggle switch connected to its context [40].(C) A narrow pulse in f 1 (input perturbation in x 1 , depicted in green) causes the isolated toggle to switch between the two stable equilibria.(D) When connected to the context, the same pulse is insufficient to yield a switch.(E) With a wider pulse, the switching is restored (however, dynamics are slower compared to the isolated case).Simulation parameters: both x 1 and x 2 bind as dimers, g 1 ~g2 ~10nM, k 1 ~k2 ~1nM 2 , d 1 ~d2 ~1hr {1 , p 1,0 ~p2,0 ~10hr {1 , g g 1 ~ g g 2 ~ g g 3 ~5nM, and the height of the input perturbation pulse is f 1 ~10nMhr {1 .doi:10.1371/journal.pcbi.1003486.g007 Remark.The above presented trade-off between robustness to interconnection and dynamic performance can be observed also in electrical systems.To illustrate this, consider the electrical circuit in Figure 9A consisting of the series interconnection of a voltage source f , a resistor R and a capacitor C, in which the output voltage is w.The speed of the circuit can be characterized by its time constant t~RC: the greater t, the slower the response.Upon interconnection with its context, represented by the capacitor C C, the time constant of the system changes to t~R(Cz C C), while the steady state remains the same.Note that the percentage change in t decreases with C, making the module more robust to interconnection, at the expense of slower response when isolated.
To further generalize the analogy between electrical systems and gene transcription networks [43], consider the electrical circuits in Figure 9B.When the module is not connected to its context, we have w~f and w w~ f f , which changes to upon interconnection.This relationship is conceptually analogous to (12).That is, the module is robust to interconnection with its context if Z is small compared to Z Z, whereas the genetic module is robust to interconnection with its context if S S and M M are ''small'' compared to R. Therefore, R is conceptually analogous to 1=Z (output admittance), whereas S S and M M play a role similar to 1= Z Z (input admittance).

Discussion
In this paper, we have focused on retroactivity, one source of context-dependence, and demonstrated that the internal, scaling, and mixing retroactivity provide missing knowledge that captures loading effects due to intramodular and intermodular connections.The internal retroactivity quantifies the effect of intramodular load, applied by nodes within a module onto each other because of binding to promoter sites within the module.Given a module of interest, the effects of intramodular loading on the module's dynamics are captured by equations ( 5)-( 6), in which one needs to As a result, the module becomes slower as R increases.(C) The percentage increase in the response time of the module decreases with g 1 , that is, internal retroactivity R increases the robustness to interconnection.Simulation parameters: d 1 ~1hr {1 , k 1 ~10nM and p 1,0 is changed such that x 1 ~50nM at the steady state.The context contains l nodes each with DNA concentration g g i ~1nM for i~1,2, . . .,l (low load: l~10; medium load: l~20; high load: l~50).The response time is calculated as the time required to reach 50% of the steady state value.doi:10.1371/journal.pcbi.1003486.g008According to the fundamental theorem by Thevenin [48], any linear electrical network can be equivalently represented by a series interconnection of a voltage source and an impedance.As a result, a generic module consists of the series interconnection of a voltage source f and an impedance Z, and similarly, any context can be represented with the series interconnection of f f and Z Z. doi:10.1371/journal.pcbi.1003486.g009 replace the specific expressions of the Hill functions H i (p i ) and node retroactivities R i (p i ) provided in Figure 3, and the binary matrices V i encoding the network topology.The scaling retroactivity accounts for the intermodular loading that the context applies on a module, due to having some TFs of the module bound to promoter sites belonging to the context.The mixing retroactivity couples the dynamics of the module and that of the context upon interconnection, and it is non-zero when TFs from different modules bind non-independently at promoter sites.The effects of intermodular loading are captured by equations ( 10)- (13).To obtain this description, it is sufficient to consider the Hill functions H i (p i ) and node retroactivities R i (p i ) provided in Figure 3, together with the binary matrices V i , D i and U representing the network topology.In general, the effects of the retroactivity matrices tend to increase with increased DNA copy number and/or decreased dissociation constants.
We have illustrated that accounting for retroactivity reveals surprising dynamical properties of modules and, at the same time, can aid design.For example, negative autoregulation, depending on the gene copy number and production rates, can slow down the response of a system instead of speeding it up.A gene can respond to a perturbation applied to a different gene even in the absence of a regulatory path between the two genes.We have shown that this can occur when a group of TFs co-regulate common targets and these common targets are found in abundance.This type of motif, referred to as the dense overlapping regulon, is highly frequent in natural regulatory networks [1].As a result, system identification techniques based on perturbation analysis [23] could erroneously identify non-existent regulatory linkages if retroactivity is not accounted for in the corresponding models.An activator-repressor clock on low copy DNA plasmid displays sustained oscillations when internal retroactivity is neglected, while oscillations are quenched once internal retroactivity is accounted for.However, by carefully adjusting the module's internal retroactivity through the addition of DNA load for the repressor, we can restore oscillations.A genetic toggle switch that can be flipped by a transient external stimulation requires a substantially longer stimulation to be flipped once it is connected to just few downstream targets.These facts are relevant, in particular, when designing synthetic biology circuits and multi-module systems.
Similar to synthetic systems, natural systems are subject to retroactivity.For example, clocks responsible for circadian rhythms have a large number of downstream targets [44], [45], which, in turn, may apply substantial load.This load can affect the amplitude and frequency of oscillations of the clock as well as the stability of the corresponding limit cycle.This suggests that natural systems may have evolved to use retroactivity in advantageous ways such as using it to properly tune the dynamic behavior of a module without changing the module's components.This hypothesis is further supported by the fact that there are a large number of TF binding sites on the chromosome that do not have a regulatory function [46], [47].These sites have an impact on the temporal response of TFs, and could therefore be exploited by nature to further control the dynamics of gene regulation.More generally, retroactivity provides means for information to travel from downstream targets to upstream regulators, therefore establishing indirect connections.In highly interconnected topologies, this information transfer can result in previously unknown ways of realizing sophisticated functions.One such example that we have provided is the adaptation function that topologies such as the dense overlapping regulon can realize by virtue of having nodes co-regulate multiple downstream targets.
Based on the three retroactivity matrices, we provided a metric of robustness to interconnection, quantifying the percent change between the dynamics of a module in isolation and once connected to other modules.This metric is an explicit function of measurable parameters and becomes smaller when a module's internal retroactivity is large compared to the scaling retroactivity of the modules it connects to.This interplay may help uncover trade-offs in natural systems, providing a new angle for understanding natural principles of network organization.From an engineering perspective, we have provided quantitative design tools that can be employed in synthetic biology to appropriately match the internal and scaling retroactivity of connected circuits to preserve the circuits' behavior upon interconnection with different contexts.Our metric of robustness to interconnection further allows to evaluate the extent of modularity of a dynamical module, possibly enabling the discovery of previously unknown core processes.Our metric could be employed by currently available methods for partitioning networks into modules.Specifically, to evaluate connectivity, these methods rely on several metrics, for instance, edge betweenness [35], its extension to directed graphs with nonuniform weights [36], round trip distance [25] or retroactivity [37].The metric of robustness to interconnection that we have introduced can enhance these methods by providing a way to evaluate modules on the basis of their functional robustness to interconnection in addition to distinguishing them at the connectivity level.
The framework that we have proposed carries substantial conceptual analogies with the electrical circuit theory established by Thevenin [48], which has been used for more than a hundred years to analyze and to design electrical networks.Within this theory, each circuit has an equivalent input and output impedance (conceptually analogous to the scaling/mixing and internal retroactivity, respectively), and an equivalent energy source (playing a role similar to the isolated module dynamics).This theory has been instrumental for answering key questions in the analysis and design of electrical networks including, for example, how the output of a circuit changes after it is interconnected in a network; how to design circuits to maximize the power transfer upon connection (impedance matching); and how to design circuits whose input/output response is unaffected by loads.We believe that the framework proposed in this paper can be used in a similar way for the analysis and design of gene regulatory networks.
Although our framework can be applied to a general gene transcription network, there are a number of aspects that it does not currently capture.These include post-translational protein modifications, such as phosphorylation, which are present in many regulatory networks and may potentially affect retroactivity.Including these will require to extend our framework to mixed gene transcription and signaling network models, leading to systems with multiple time scales.Furthermore, the transcription and translation processes use shared resources such as RNA polymerase and ribosomes, which may create couplings among unconnected circuits [17].The dynamics of shared resources has not been included in our modeling framework and will be the focus of our future work.

Detailed Description of the System Model
The production of TF x i is regulated by its parents p i,1 ,p i,2 , . ..: they bind to the promoter of x i , and form complexes c i,1 , c i,2 , . . .with the promoter.Each of these complexes, in turn, produce x i with a different rate, where we use a one-step production process encapsulating both transcription and translation [3].As a result, the reactions we consider for node x i are modeling the following physical phenomena.We denote by d i protein decay, whereas f i t ð Þ represents the production rate that may be due to external inputs or perturbations (inducer, noise or disturbance).The second reversible reaction in (17) describes the binding of parent p i,l with multimerization factor m i,l to promoter complex c i,j forming complex c i,k , where a i,j,k and b i,k,j are the association and dissociation rate constants, respectively.Furthermore, each promoter complex c i,j will contribute to the production of x i through the production rate constant p i,j , modeled by the third reaction in (17).This production rate constant is a lumped parameter that incorporates features such as the RBS strength and the promoter strength.Finally, we assume that the total concentration of the promoter, denoted by g i , for each transcription component is conserved, so that g i ~PC i j~0 c i,j , where C i is the number of possible complexes formed with the promoter of x i .This concentration is proportional to the concentration of copies of the promoter, which can be controlled, for example, by changing plasmid copy numbers in synthetic systems.
The reaction flux vector v contains all the reactions in the system, that is, binding/unbinding and protein production/decay.Given that binding/unbinding reactions occur on a much faster time-scale than protein production/decay [1], we partition v into r Ã and r, where r Ã contains the slow processes, whereas r is composed of the fast reactions, that is, r Ã ~. .

Biochemical Parameters
Since the production of a typical protein takes approximately 5 minutes [1], and a few dozen mRNAs can be transcribed from the same gene simultaneously by [49], and similarly, a few dozen proteins can be translated from the same mRNA at the same time by [50], the effective production rate of protein from a gene can be as high as p&10000hr {1 .This value can be arbitrarily decreased, for instance, by decreasing the RBS strength in synthetic circuits.The cell volume of Escherichia coli is typically between 0:34{1:32mm 3 by [51], so that 1 molecule/cell corresponds to approximately 1{5nM concentration.By [52], a typical value of the dissociation constant of bacterial promoters is k~1nM, whereas [22] suggests k~10nM, and experimentally obtained values are provided in [53].One of the most widely used high copy number vectors is the pUC plasmid [54], which can have hundreds of copies per cell [55].A frequently used medium copy number plasmid is p15A with a few dozen copies per cell [56], whereas pSC101 is regarded as a low copy number plasmid with only a few copies per cell [56].Finally, since the lifetime of a protein is on the order of a cell-cycle [1], we have d~0:3{1:2hr {1 [50].The typical range of macroscopic parameters in Escherichia coli is summarized in Table 1.
If we had not neglected mRNA dynamics, there would be three different time scales in the system.Binding and unbinding reactions occur on the time scale of seconds (or even subseconds) [1], representing the fastest time scale.The intermediate time scale is that of mRNA dynamics, as the average lifetime of mRNA is on the time scale of minutes [57], [58], [59].Finally, the dynamics of proteins evolve on the slowest time scale (hours).As we are interested in describing the dynamics of the system on the time scale of gene expression, the concentration of promoter complexes and mRNA transcripts can be both approximated with their quasisteady state values, leading to the models we have proposed in the paper.However, we would like to point out that including mRNA dynamics would not change anything substantial in the results and it would simply add N more ODEs to the ODE model of an N-node module without any effects on the retroactivity matrices (shown in [8] considering a specific example).
Definition of H i (p i ) and R i (p i ) where r i (p i ,c i ) denotes the reaction flux vector corresponding to reversible binding reactions with the promoter of x i .Let c i ~ci (p i ) denote the vector of concentrations of complexes with the promoter of x i at the quasi-steady state, obtained by setting 0~A i r i (p i ,c i ).
We first define H i (p i ) as follows: Table 1.Typical range of macroscopic parameters in Escherichia coli.

Figure 2 .
Figure2.The context (downstream system) affects the behavior of the module (upstream system).(A) The module in isolation.(B) The module in isolation displays sustained oscillations.(C) The module connected to its context.(D) Upon interconnection with its context, the dynamics of the module change due to the retroactivity s from its context, since some of the molecules of x 1 are involved in binding reactions at node x x 2 .As a result, those molecules are not available for reactions in the module, and the output of the module is severely changed.For details on the system and parameters, see Supporting Text S1. doi:10.1371/journal.pcbi.1003486.g002

Figure 3 .
Figure3.Hill function and retroactivity of node x i for the most common binding types.If node x i has no parents, its node retroactivity is not defined.In the single parent case, node x i has one parent, y binding as an n-multimer with dissociation constant k y .In the case of independent, competitive and cooperative binding, node x i has two parents, y and z, binding as multimers with multimerization factors n and m, respectively, together with dissociation constants k y and k z , respectively.The total concentration of the promoter of x i is denoted by g i .The production rates p i,0 , p i,1 , p i,2 and p i,3 correspond to the promoter complexes without parents, with y only, with z only, and with both y and z, respectively.For details, see Supporting Text S3. doi:10.1371/journal.pcbi.1003486.g003

Figure 4 .
Figure 4. Negative autoregulation can make the temporal response slower.Time response at a steady state fixed at x 1 ~50nM.The red and blue plots denote the cases with and without negative autoregulation, respectively, whereas the green plot represents the case of negative autoregulation neglecting retroactivity (R~0 in(8)).Simulation parameters are d 1 ~1hr {1 , k 1 ~10nM, together with p 1,0 ~20hr {1 , p 1,0 ~10hr {1 , p 1,0 ~1hr {1 for A, B and C, respectively.To carry out a meaningful comparison between the unregulated and regulated systems, we compare the response time of systems with the same steady state.To do so, we pick the same value of g 1 in the case of the regulated systems (g 1 ~15nM, g 1 ~30nM, g 1 ~300nM for A, B and C, respectively), but a different one for for the unregulated system (g 1 ~2:5nM, g 1 ~5nM, g 1 ~50nM for A, B and C, respectively), such that the steady states match (see Methods for parameter ranges).Decreasing p 1,0 (lower production rate constant) while increasing g 1 (higher DNA copy number) results in slower response, as internal retroactivity increases.doi:10.1371/journal.pcbi.1003486.g004

Figure 5 .
Figure 5. Nodes can become coupled via common downstream targets.(A) Node x 3 has two parents: x 1 and x 2 , without a regulatory path between them.(B) Perturbation f 1 applied to x 1 .(C) In the case of competitive binding, increasing the concentration of free x 1 yields more of x 1 bound to the promoter of x 3 , forcing some of the molecules of x 2 to unbind, thus increasing the free concentration x 2 .Consequently, x 1 acts as if it were an activator of x 2 .(D) By contrast, in the case of cooperative binding, when the binding of x 1 must precede that of x 2 , pulses in x 1 yield pulses of the opposite sign in x 2 .Consequently, x 1 acts as if it were a repressor of x 2 .Simulation parameters are: g 1 ~g2 ~10nM, g 3 ~20nM, d 1 ~d2 ~1hr {1 , p 1,0 ~0, p 2,0 ~10hr {1 , k 1 ~k2 ~1nM 4 , and both x 1 and x 2 bind as tetramers.doi:10.1371/journal.pcbi.1003486.g005

Figure 7 .
Figure 7. Effects of the context on the switching characteristics of the toggle switch.(A) The toggle switch in isolation.(B)The toggle switch connected to its context[40].(C) A narrow pulse in f 1 (input perturbation in x 1 , depicted in green) causes the isolated toggle to switch between the two stable equilibria.(D) When connected to the context, the same pulse is insufficient to yield a switch.(E) With a wider pulse, the switching is restored (however, dynamics are slower compared to the isolated case).Simulation parameters: both x 1 and x 2 bind as dimers, g 1 ~g2 ~10nM, k 1 ~k2 ~1nM 2 , d 1 ~d2 ~1hr {1 , p 1,0 ~p2,0 ~10hr {1 , g g 1 ~ g g 2 ~ g g 3 ~5nM, and the height of the input perturbation pulse is f 1 ~10nMhr {1 .doi:10.1371/journal.pcbi.1003486.g007

Figure 8 .
Figure 8. Internal retroactivity a module more to at the expense of speed.(A)The module consists of a single negatively autoregulated node, whereas the context comprises l nodes repressed by the TF in the module.(B) The internal retroactivity R of the module increases with the DNA copy number g 1 of x 1 .As a result, the module becomes slower as R increases.(C) The percentage increase in the response time of the module decreases with g 1 , that is, internal retroactivity R increases the robustness to interconnection.Simulation parameters: d 1 ~1hr {1 , k 1 ~10nM and p 1,0 is changed such that x 1 ~50nM at the steady state.The context contains l nodes each with DNA concentration g g i ~1nM for i~1,2, . . .,l (low load: l~10; medium load: l~20; high load: l~50).The response time is calculated as the time required to reach 50% of the steady state value.doi:10.1371/journal.pcbi.1003486.g008

Figure 9 .
Figure 9. Analogy with electrical systems.(A) The module consists of the series interconnection of a voltage source f , a resistor R and a capacitor C. The speed of the module can be characterized by the time constant t~RC, which increases upon interconnection with the context.The greater C, the slower the module in isolation, but the smaller the percentage change in its speed upon interconnection.(B)According to the fundamental theorem by Thevenin[48], any linear electrical network can be equivalently represented by a series interconnection of a voltage source and an impedance.As a result, a generic module consists of the series interconnection of a voltage source f and an impedance Z, and similarly, any context can be represented with the series interconnection of f f and Z Z. doi:10.1371/journal.pcbi.1003486.g009

First, note
that A in (2) has a block diagonal structure, yielding _