Retroactivity to the input in complex gene transcription networks

Synthetic biology is a bottom-up engineering discipline: biological modules are systematically designed with predefined behavior and then combined to build up larger circuits. Although the modules produce the desired behavior in isolation, they fail to operate properly when they are connected due to retroactivity, an effect which extends the notion of impedance to biomolecular systems. Despite playing a central role, retroactivity is not yet characterized in complex gene transcription networks. In this paper, we mathematically describe and quantify this effect. This result is obtained by applying singular perturbation on the finite time interval. We identify the biomolecular counterpart of impedance and introduce the effective retroactivity to the input of a gene. Furthermore, we provide a theorem describing how modules affect each other when connected. We restore modular composition of synthetic circuits by extending the characterization of modules with internal and input retroactivities. We illustrate the implications of the results by investigating crosstalk in a simple genetic system.


I. INTRODUCTION
Modularity is a central concept in every engineering discipline. In lack of it, building large, complex systems by composing smaller, simple pieces together is infeasible. Digital electronics prominently illustrates this: with functional modules such as clocks, memories and arithmetic units one can design large systems by connecting these basic building blocks. What makes this approach powerful is the fact that modules behave the same as if they were in isolation.
In case of analog electronics, modules typically affect each other. However, the behavior of connected modules can be easily described by considering their models in isolation. The fundamental theorem by Thevenin [1] makes it possible to substitute electrical circuits between any two terminals by an equivalent circuit consisting of the series connection of a single voltage source and impedance. When connecting modules, one can consider the equivalent models describing the modules' behavior in isolation. This result heavily relies on the fact that the impedance of an electrical component remains unchanged when connected to other components.
Synthetic biology is closely related to analog electronics. The basic building blocks are transcription components producing a single transcription factor (TF) as output and taking a few TFs as input. Instead of wires, transcription components are connected via binding reactions: input TFs bind to the promoter region and as a result of transcription and translation processes the output TF is produced. Given the close relationship between synthetic biology and analog electronics, it is natural that an impedance-like effect is observable when connecting biomolecular components together: this effect is called retroactivity [2]. Retroactivity arises whenever two molecules bind together describing the fact that these molecules become unavailable for other reactions. A key feature of retroactivity is that it enables a downstream system to affect the behavior of the upstream one [3], [4]. In spite of its central role, retroactivity is not yet characterized for complex gene transcription networks.
Therefore, in this paper we characterize retroactivity in gene transcription networks with arbitrary topology. We define the effective retroactivity to the input of a transcription component and we argue that it can be interpreted as the biomolecular analog of impedance. We introduce the internal retroactivity of a module capturing the effect of intramodular connections. This is followed by our main result: a theorem for complex gene transcription networks describing how the dynamics of modules change upon interconnection. We introduce the effective retroactivity to the input of a module, a quantity similar to input impedance. We show how the dynamics of interconnected modules can be determined considering (i) their dynamics in isolation, (ii) their internal retroactivity and (iii) their retroactivity to the input. We therefore recover a modular approach to understand the dynamics of complex systems by augmenting their description with internal and input retroactivities. For the most common binding types (independent, cooperative and competitive) we provide the explicit expression of the effective retroactivity to the input of a transcription component. This means that having a transcription network where the binding reactions are of these basic types, one can compute the internal and input retroactivity of a module just as easily as in case of electrical circuits. In order to show the power of the framework, we investigate crosstalk between modules.
Our work is complementary to those partitioning large transcription networks into modules by minimizing retroactivity ( [5], [6], [7] and [8]). Here, we analytically characterize and quantify retroactivity using singular perturbation theory. Singular perturbation has been used before as a powerful model reduction arpproach for gene network models [9]. The notion of retroactivity connects with the idea of fan-out introduced in [10]. Our approach is based on the tools of dynamical systems analysis, hence it connects with other disciplines of biochemical systems analysis, such as metabolic control analysis [11], [12] and metabolic supply and demand analysis [13]. However, whereas these methods are primarily interested in the steady state and near-equilibrium behavior of a system, we focus on the dynamics of modules and biomolecular circuits.

II. SYSTEM MODEL
We view a gene transcription network as a hierarchical structure with three levels: nodes represent transcription components, a group of connected nodes forms a module, whereas a system consists of several modules. Throughout the paper, species are denoted by capital letters, whereas the corresponding lower case letter stands for their concentration, e.g., the concentration of free TF X 2 is x 2 . Moreover, let us use the superscript for referring to modules, that is, X M 2 belongs to module M .

A. First Hierarchical Level: Nodes
Transcription networks are usually viewed as the input/output interconnection of fundamental building blocks called transcription components (nodes). A transcription component (Fig. 1a) takes a number of TFs as input forming complexes with the promoter sites through reversible reactions to produce a single TF as output, through the process of gene expression [14].
Denote by Z p,k (k = 1, 2, . . . , ζ p ) the TFs regulating X p , and call them the parents of node p (Fig. 1b). Moreover, define Z p as the set of parents of node p. The concentration of Z p,k is denoted by z p,k , whereas z p stands for the concentration vector of parents of node p, that is, z p z p,1 z p,2 . . . z p,ζ p T .
Furthermore, denote by C p,i (i = 0, 1, 2, . . . , χ p ) the possible complexes formed by the promoter of X p with some Z p,k ∈ Z p such that C p,0 is the free promoter. Moreover, let C p be the set of possible complexes corresponding to node p and define c p as the concentration vector of complexes at node p except for the free promoter, i.e., c p c p,1 c p,2 . . . c p,χ p T .
The reactions we consider for node p are as follows. We denote by ν p and δ p the external input and protein decay, respectively, that is, Reversible binding reactions are characterized by rate constants α p,i , j and α p, j ,i Node p with its parents Z p Fig. 1: (a) The promoter contains the regulatory sites where input TFs (X q , X r , . . . , X s ) can bind forming complexes C p,i (i = 0, 1, . . . , χ p ). The coding region encodes the genetic information required for the expression of the output TF (X p ). (b) There is a directed edge from node q to p if X q is an input to the transcription component producing X p , and we refer to X q as the parent of node p.
we denote by π p,i the rate constants describing protein production: C p,i Finally, we assume that the total concentration of promoter (η p ) is conserved: Define the set Ω p as follows: i , j , k ∈ Ω p if TF Z p,k can bind to complex C p,i forming complex C p, j , i.e., For instance, consider node p having two parents Z p,1 and Z p,2 binding cooperatively, that is, first Z p,1 has to bind to the free promoter C p,0 forming complex C p,1 , and only after that can Z p,2 bind resulting in complex C p,2 . In this case Considering the reactions for node p, one can writė In addition, it follows from (1) that Consequently, we can disregard one of the equationṡ c p,i = ϕ p,i c p , z p , for instance the one standing for the free promoter (i = 0). Introducing ϕ p c p , z p ϕ p,1 c p , z p ϕ p,2 c p , z p . . . ϕ p,χ p c p , z p T , the evolution of complexes at node p is described bẏ If X p is not taken as input by any nodes (including node p itself ), that is, X p does not take part in any binding reaction, then its dynamics can be described bẏ Assume now that X p is taken as input to other nodes, that is, X p takes part in binding reactions. Having a total of n nodes, define c c T concentration vector of all complexes (except for the free promoters). Denote by b p (c) the concentration of bound X p , or equivalently, the concentration of complexes having Note thatḃ p (c) represents the rate of change of bound X p . Consequently, the rate of change of free X p due to binding reactions is −ḃ p (c). Combining this with the reactions considered in the system, we obtainẋ p = f p x p , c with Definition 2. The connected dynamics of node p are defined asċ p = ϕ p c p , z p andẋ p = f p x p , c .

B. Second Hierarchical Level: Modules
Modules are considered to be a group of connected nodes with some functionality.
. . , X n } is the set of TFs in the module; • C = {C 1 ,C 2 , . . . ,C n } is the set of complexes in the module, where C p = C p,0 ,C p,1 , . . . ,C p,χ p is the set of complexes formed at node p; where U i is a TF from a different module; • P = {P 1 , P 2 , . . . , P n } is the set of parameters describing the reactions of the module, where P p is the set of parameters associated with node p.
Since a module is an ensemble of connected nodes, the ODE model of a module is simplyċ p = ϕ p c p , z p anḋ . . x n T , the concentration vector of free TFs in the module and f (x, c) Note that in (4) the argument on the right hand side is c p and z p for p = 1, 2, . . . , n, whereas on the left hand side it is x, c and u. This is because ∪ n p=1 Z p ⊆ X ∪U , i.e., parents in the module are either nodes in the module or inputs. The above definition describes the case when TFs of the module are not taken as input to any other module, that is, the module is in isolation.
Furthermore, the module's dynamics without considering the loading effect of intramodular binding reactions are given byċ = ϕ (x, c, u) andẋ = g (x, c) with g (x, c) being the column vector of g p x p , c p for p = 1, 2, . . . , n. Example 1. Consider the system in Fig. 2, which will serve as a running example. Given that X N 1 has two parents, we choose Z N 1,1 = X M 1 and Z N 1,2 = X N 1 . There are four complexes associated with node 1 in module N : the free promoter at this node is denoted by C N 1,0 , whereas C N 1,1 and C N 1,2 stand for the complexes of promoter with X M 1 and X N 1 , respectively. Finally, C N 1,3 denotes the complex of promoter with both TFs bound.

C. Third Hierarchical Level: Systems
Definition 5. We say modules M and N are composable if X M ∩ X N = , that is, the modules do not share nodes.

III. EFFECT OF INTERCONNECTIONS
Our first question relates to connecting nodes: what is the relation between the isolated and connected dynamics of a node? The second question focuses on connecting modules: how do the dynamics of composable modules change upon interconnection? For simpler notation, we only use the superscript when we need to distinguish modules, i.e., when there are multiple modules in focus.

A. Effective Retroactivity to the Input of a Node
Consider n interconnected nodes, that is, a module. Define the parent matrix of node p as Ψ p = Ψ p ζ p ×χ p , where Ψ p k, j = 1 if there exists i such that i , j , k ∈ Ω p , otherwise it is zero. Since Ψ p k, j = 1 means that complex C p, j has TF Z p,k bound, we can calculate the concentration of bound parents at each node by defining with w p c p = w p,1 c p w p,2 c p . . . w p,ζ p c p T .
Note that w p,k c p denotes the total concentration of complexes at node p having TF Z p,k bound, or equivalently, the total concentration of bound Z p,k at node p (p = 1, 2, . . . , n and k = 1, 2, . . . , ζ p ). Let γ p z p = γ p,1 z p γ p,2 z p . . . γ p,χ p z p T be the solution of 0 = ϕ p γ p z p , z p , the concentration of complexes at node p whenċ p = 0. By (1) we can define γ p,0 z p η p − χ p i =1 γ p,i z p , the concentration of C p,0 (free promoter) whenċ p = 0.
Denote by Φ the set of nodes in the module having parents, that is, Φ = p Z p = .
Definition 7. Define the effective retroactivity to the input of node p as In other words, R p z p denotes the sensitivity of the total concentration of bound parents to the concentration of free parents at node p when c p = γ p z p , that is, wheṅ c p = 0. Furthermore, by (6)-(7) we obtain Since γ p z p only depends on parameters of node p, R p z p also depends only on the parameters of node p. Therefore, R p z p is the property of the node and it is independent of network topology, that is, it does not change upon interconnection. Furthermore, one can verify that R p z p is the generalization of retroactivity introduced in [2] for combinatorial regulation.  (9) and considering the regulation patterns in Example 1: Making later computations simpler, let us write where a, b, . . . , f are implicitly defined in (9)- (12).

B. Effect of Intramodular Connections
Here, we show that the isolated dynamics of module M given in Definition 4 can be well approximated by considering the isolated dynamics of nodes in M and the effective retroactivity to the input of node p for p ∈ Φ M .
Consider node p ∈ Φ M (a node in M having parents) and some module N . Define the n N ×ζ M p transformation matrix p,N provides us with a mapping between parents of node p ∈ Φ M and nodes in N .
Based on (7), every row and column in R M p z M p correspond to a parent Z M p,k ∈ Z M p of node p ∈ Φ M . Every parent Z M p,k ∈ Z M p is either a node in N or not. With the restriction operator we select the rows and columns of R M p z M p corresponding to parents belonging to N , and rearrange them according to the order of nodes in N .
from Example 3 with (13). Its first row and column belong to Z N 1,1 = X M 1 ∈ X M (parent from M ), whereas the second row and column are associated with Z N 1,2 = X N 1 ∈ X N (parent from N ). Therefore, when considering the restriction of R N 1 z N 1 to M , we select b. Since X M 1 is the only node in M , we have R N 1 z N 1 M = b. In case of R N 1 z N 1 N , we select e and since it belongs to parent X N 1 which is the first node in N , we obtain R N 1 z N 1 N = e 0 0 0 . Similarly, when considering the restriction of R N 1 z N 1 to M N , we select the whole matrix and rearrange its rows and columns according to the the order of nodes in M N and obtain

Definition 8. The internal retroactivity of module M is
where Φ is the set of nodes in M having parents.
According to [14], the binding reactions are much faster than protein production and decay. Therefore, picking protein decay rate δ r = 0 and dissociation rate α s,k,l = 0 is a dimensionless small parameter, that is, ε 1.
provided that the matrix ∂ϕ ∂c c=γ (x,u) is Hurwitz.
where ε is the small parameter from (16). Consequently, (4)- (5) becomė which is in the standard singular perturbation form, where ξ is the slow variable, whereas c is the fast variable. By setting ε = 0 in (19), we obtain the slow manifold [15] on which the dynamics of the system are governed by the slow variable dynamics. It can be shown that in this case c = γ (x, u) andξ = g x, γ (x, u) . Furthermore, = R (x, u), thusẋ =f (x, u) describes the dynamics of (4)-(5) on the slow manifold. Since we assume that ∂ϕ ∂c c=γ (x,u) is Hurwitz, the slow manifold is locally exponentially stable, hence the dynamics restricted to the slow manifold are a good approximation (Theorem 11.4 in [15]), which completes the proof. ■ Theorem 1 states thatx (t ) approximates x (t ) well if ε 1, thus we refer toẋ =f (x, u) as the reduced order model of module M in isolation.
Looking at (17), R (x, u) relates the dynamics of the connected and isolated nodes in M . In other words, R (x, u) captures the retroactive effects due to intramodular binding reactions, hence the notion internal retroactivity.
Recalling that R p z p does not change upon interconnection, and that it captures the loading effect from downstream nodes, it can be interpreted as the biomolecular analog of impedance.
is white noise with unit power, ν N 2 (t ) ≡ 0 and all the production rate constants are zero.
In case of cooperative binding we have c, d = 0 by (11). Therefore, the off-diagonal terms in (I +G∆) −1 are nonzero. It follows that we have to consider two effects. First, the isolated module dynamics of x M 1 and x N 1 are scaled by the first two diagonal terms when connected. Second, an additive crosstalk occurs between x M 1 and x N 1 through the off-diagonal terms. One might think that the latter effect from M to N is because M is upstream whereas N is downstream, but this is clearly incorrect: in case of independent binding we do not have this phenomenon. The crosstalk here is purely due to the nonindependent (cooperative) binding and its extent is determined by the magnitude of the off-diagonal terms.
Simulation results for cooperative binding (Fig. 3) confirm that the isolated behavior of the downstream module N distorts the periodic output signal of the upstream system M . This will always occur whenever c = 0, that is, the binding of X M 1 is not independent of the binding of X N 1 . The larger c(1+e+ f ) (1+a+b)(1+e+ f )−cd , the greater this effect. Taking (11), one can see that c ≈ 0, for instance, if η N 1 ≈ 0 or x M 1 ≈ 0, that is, if the total concentration of the downstream system or the output signal of the upstream system is small. Furthermore, we have c ≈ 0 if the dissociation constant k M is large. All these conditions represent that the connection between M and N is weak, thus the crosstalk from N to M is small.

V. CONCLUSION AND FUTURE WORK
In this paper, we applied singular perturbation theory to study retroactivity and modularity in complex gene transcription networks.
First, we introduced R p z p , the effective retroactivity to the input of node p and argued that it can be interpreted as input impedance. It only depends on parameters associated with the node, that is, it remains unchanged when the transcription component is part of a larger network. Furthermore, it describes the loading effect when a downstream component is connected. In addition to providing a formula for calculating this key quantity, we presented the expression of R p z p for the most common regulation types: independent, cooperative and competitive.
Second, we defined the internal retroactivity of a module capturing the retroactive effects due to intramodular connections. Moreover, we introduced a module's effective retroactivity to the input to another module describing the load presented by intermodular binding reactions when connecting two modules.
Finally, we presented a theorem for complex gene transcription networks analog to Thevenin's. It allows us to determine the behavior of connected modules by considering (i) their model in isolation, (ii) their internal retroactivity and (iii) their effective retroactivity to the input to each other.
Although the current framework is capable of modeling the most relevant processes, such as protein production and decay, as well as binding and unbinding reactions, we will extend our approach by including mRNA dynamics and dimerization. In addition, we propose to investigate the effect of retroactivity for complex systems from a qualitative point of view.