Robust Signal Processing in Living Cells

Cellular signaling networks have evolved an astonishing ability to function reliably and with high fidelity in uncertain environments. A crucial prerequisite for the high precision exhibited by many signaling circuits is their ability to keep the concentrations of active signaling compounds within tightly defined bounds, despite strong stochastic fluctuations in copy numbers and other detrimental influences. Based on a simple mathematical formalism, we identify topological organizing principles that facilitate such robust control of intracellular concentrations in the face of multifarious perturbations. Our framework allows us to judge whether a multiple-input-multiple-output reaction network is robust against large perturbations of network parameters and enables the predictive design of perfectly robust synthetic network architectures. Utilizing the Escherichia coli chemotaxis pathway as a hallmark example, we provide experimental evidence that our framework indeed allows us to unravel the topological organization of robust signaling. We demonstrate that the specific organization of the pathway allows the system to maintain global concentration robustness of the diffusible response regulator CheY with respect to several dominant perturbations. Our framework provides a counterpoint to the hypothesis that cellular function relies on an extensive machinery to fine-tune or control intracellular parameters. Rather, we suggest that for a large class of perturbations, there exists an appropriate topology that renders the network output invariant to the respective perturbations.


I. INTRODUCTION
We aim at a mathematical formalism that allows to judge whether a multiple-inputmultiple-output reaction network is robust against large perturbations of network parameters. In addition to identify design principles for robust signal processing, the formalism should indicate the necessary network modifications and extensions to arrive at a robust network output for specific perturbations, when starting from a non-robust network. In particular, our formalism builds upon the structural properties of a (bio)chemical network, as it is the network architecture and not 'fine-tuning' of parameters that allows for the compensation of large perturbations. Our results show that the robustness of a signaling network can be judged by inspection of a linear vector space: We demonstrate that for each biochemical network, there exists a linear vector space, such that any perturbation (expressed as a vector of partial logarithmic derivatives) that is confined to this vector space leaves the output of the network, as defined by a set of stationary concentrations of designated output variables, invariant. One of our main achievement is to identify which part of this vector space is independent of kinetic parameters (corresponding to global concentration robustness) and which part of this vector space is dependent on kinetic parameters (corresponding to local concentration robustness and requiring fine-tuning of parameters). The former is denoted as the invariant perturbation space of the network and can be algorithmically constructed for any signaling network.
Our framework provides a counterpoint to the hypothesis that cellular function relies on an extensive machinery to fine-tune or control intracellular parameters. Rather, our framework suggests that there exists an appropriate topology that renders the network output insusceptible to a given class of perturbations. Our framework draws upon and extends concepts of (metabolic) control theory. The reader familiar with nonlinear control theory will discover some parallels to the concept of disturbance decoupling. Within this Supplementary Information, we first outline our framework in more detail, point out pitfalls, and provide additional examples. A mathematical rigorous treatment is given in Section IV.

A. Biochemical networks without conserved moieties
We first illustrate our concept using the special case of a biochemical network without conserved moieties. The biochemical network is assumed to consist of m independent dynamic state variables, x = (x 1 , ..., x m ), whose temporal evolution is determined by a differential with α ij denoting the kinetic exponents that do not necessarily have to take integer values.
With respect to the state variables, we further distinguish between a set of output state variables, defined as x A , and a set of intermediate state variables, x M . As the prerequisite for robustness, we require that the output states are invariant under perturbations, that is, ∆x A = 0, where ∆x A denotes the difference in the output variables after and before a perturbation ∆p j on network parameters.
In the following, we assume the existence of a (not necessarily unique) asymptotically stable stationary state x s with N · v(x s ) = 0. We further assume that the functionality of the network is encoded in the stationary dependence of the designated output variables, x A , on a set of input parameters.
We expand the stationary form of the Eq. (1), N · v s = 0, with v s := v(x s ), to linear order in both the perturbations, ∆p j , that change network fluxes v s i (p j ) → v s i (p j + ∆p j ) and the resulting changes in the state variables ∆x, with diag(v s ) denoting a square matrix with entries v s on the diagonal and the expansion coefficients The relative perturbations and its responses are defined as (∆p We note that if the reaction fluxes follow the functional form given in Eq. (2), the expansion coefficients are given by the constant kinetic coefficients α ij .
In the absence of the condition ∆x A = 0, the expansion Eq. (3) always has a unique solution ∆x that quantifies the local linear response to an infinitesimal perturbation in parameters. The existence of the solution is guaranteed by the fact that the Jacobian of the system is of full rank and hence invertible -a condition that is extensively utilized within, for example, Metabolic Control Analysis [3,7].
However, our requirement of robustness ∆x A = 0 removes those degrees of freedom that correspond to (changes in) the output variablesx A . As a consequence, only the set of intermediate variablesx M is able to compensate the perturbations. In this case, Eq. (3) translates into the condition In general, Eq. Or, the columns of the matrix P are linearly dependent on the columns of the matrix M .
In mathematical terms, these two conditions can be summarized in the equation rank(P |M |K) = rank(M |K) .
Here, the columns of K span the right nullspace of N · diag(v s ), such that N · diag(v s ) · K = 0. The notation (M |K) denotes a concatenation of the columns of the matrices M and K. Equation (6) expresses the condition that each column vector of the matrix P must be linearly dependent on the column vectors of (M |K).
We note that, considering again the full system in Eq. (3), the concatenated matrix (A|M |K) is square and of full rank. This property again reflects the fact that there exists a local linear response in the systems variables for any local perturbation in the rate equations -a consequence of the Jacobian being invertible. Conversely, the concatenation (M |K) is not of full rank, hence Eq. (6) cannot be fulfilled for arbitrary perturbations.
However, if Eq. (6) is fulfilled, then the system indeed exhibits local concentration robustness, as defined by ∆x A = 0. Here, we emphasize that the linear system Eq. In addition to the rank condition Equation (6)  As yet, the conditions for local concentration robustness were derived with respect to infinitesimal perturbations at a particular state x s . In general, the matrices M and K will depend on the particular state at which the expansion was performed -hence Eq. (6) does not represent a sufficient condition for global robustness. In the following, as one of the major achievement of our work, we will derive the architectural requirements on the network that ensure that Eq. (6) is fulfilled independent of the particular stationary state, hence the system allows for global concentration robustness in the face of perturbations of large magnitude. Prior to this step, we briefly extend our analysis to networks that incorporate mass-conservation relationships.
Most models of biochemical networks exhibit conserved moieties that usually arise from an approximation of slowly changing components by constant quantities. In this case, the system of differential equations for the independent state variables, x = (x 1 , ..., x m ) is augmented by a set of dependent state variables x D , whose values are determined by n mass conservation equations. The full system of equations governing the time evolution of the system isẋ with expansion coefficients The relative perturbations and its responses are again defined as (∆p We emphasize that the changes in states ∆x D are entirely determined by changes in the intermediate states, ∆x M , via the differential mass conservation equation. Consequently, the dependent state variables do not represent additional degrees of freedom within the system and are not able to compensate perturbations.
Rather, the associated matrix D must be considered as indirect perturbations on the network that are induced by changes ∆x M .
In mathematical terms, we can use the differential mass conservation relationship to substitute changes in the dependent variables by changes in independent intermediate variables, and L ′ is obtained from L by removing those columns that correspond to output variables.
The condition for invariance of the output Eq. (6) can then be written as Eq. (12) is the generalized rank condition for local invariance of the output variables with respect to infinitesimal perturbations. A general definition of the matrix M is thus given In absence of conservation equations, L = 0, we obtain the identity M = M D and thus recover our previous result Eq. (6).
As observed in Eq. (13), mass conservation relationships induce additional elements (dependencies) in the matrix of partial derivatives -a consequence of the substitution of the dependent variables within the kinetic rate equations.
We illustrate this point with a simple example. Consider a canonical two-component signal transduction network in bacteria, e.g. the EnvZ/OmpR system, where the histidine kinase with total concentration Z T gets autophosphorylated and transfers the phospho-group to the response regulator. The response regulator, with total concentration R T , is in turn dephosphorylated proportional to the phosphatase activity of the histidine kinase, which is The system of differential equations for the independent variables is given as The reaction network further satisfies the conservation equations using the assignments x T = (Z T , R T ), x = (Z p , R p ) and x D = (Z, R), and assuming R p as the output variable of the pathway.
Utilizing the definitions given above, we then obtain hence the matrix M is given as, using the definition α := −Z p /Z = −Z p /(Z T − Z p ). As compared to the situation without mass conservation relationships the matrix M contains additional elements, corresponding to the implicit dependencies of the dependent variables. For example consider the element In Section VII A of this Supplementary Information we will consider two-component systems and robustness against variation in total compound concentrations in more detail.

C. From local to global robustness: The invariant perturbation space
The basis of our approach is the transition from local concentration robustness to a criterion of global concentration robustness. In this respect, we require that the rank condition for local concentration robustness is fulfilled at any stationary state and irrespective of the kinetic parameters. Consequently, a perturbation that is large in magnitude may gradually alter the stationary state of the system by affecting the set of intermediate variables x M -however, the stationary concentrations of the set of designated output variables x A are not affected. We note that in the following, unless otherwise stated, we always assume that the system gives rise to a globally stable steady state for all parameters. As a consequence, we assume the Jacobian to be invertible at each point in state-space. For a more formal treatment see also Section IV.
The condition for local concentration robustness is given by Eq. (6), The task is to ascertain whether the equation is fulfilled independent of kinetic parameters.
To this end, we recall that the rank of a matrix is unchanged under elementary matrix operations (EMO), which are: (i) the exchange of any two columns (rows), (ii) the multiplication of a column (row) by the same non-zero factor, (iii) the addition of an arbitrary multiple of one column (row) to another. Utilizing a series of EMOs, we aim to remove explicit parameter dependencies from Eq. (6), thus obtaining a global structural condition for concentration robustness.
To exemplify the application of EMO we continue with the example, a canonical twocomponent system, discussed in Eq. (14). In addition to the matrix M , defined above, we construct the right nullspace K of the scaled stoichiometry, consisting of one column vector K = (β β β) T , with β := (v 1 ) −1 . We are interested in concentration robustness of the output variable R p with respect to a perturbation p that affects both conserved total concentrations, Z T = Z T (p) and R T = R T (p). In this case, the perturbation vector reads with γ := ∂ ln v 1 /∂ ln p = ∂ ln v 3 /∂ ln p, and δ := ∂ ln v 2 /∂ ln p (but see also subsequent sections for introductory examples). The rank condition is therefore given as with M and α defined above. A parameter-independent representation is obtained by EMO Obviously, the columns of (M |K) span a two dimensional plane that does not change its orientation in the three dimensional space under changes in v s , x M , and kinetic parameters that enter the equations via the expressions for α, β, γ, and δ. Thus any vector P that lies within this plane fulfills the rank condition. In the example considered here the vector P = (γ δ γ) T lies in the plane spanned by column vectors (0 1 0) T and (1 0 1) T .
Since the rank condition is fulfilled irrespective of the stationary state and kinetic parameters, the stationary network output, the variable R p is globally robust with respect to P .
However, in general not all dependencies on the stationary state of the matrix (M |K) can be removed by elementary matrix operations (EMO). We therefore define by I the space spanned by the largest possible set of parameter independent column vectors of (M |K), with (M ′ |K ′ ) a reduced, maximally parameter free representation. We call I the invariant perturbation space. The invariant perturbation space contains exclusively structural information of a reaction network. Our approach therefore allows to separate the structural from fine tuned network properties. In terms of robustness -as defined in this work -this means that any perturbations P of the reaction fluxes that lie entirely in the invariant subspace result in an invariant system output. If we define by (P ′ |M ′ |K ′ ) the matrix resulting from (P |M |K) after appropriate elementary matrix operations, then is a sufficient condition for structural robustness. For a rigorous definition see Section IV.

D. Modelling biochemical reaction networks
Our approach relies on an interpretation of the network structure in terms of the logarithmic partial derivatives of the kinetic rate equations. In particular, we utilize the fact that the logarithmic partial derivatives are -for certain kinetic functions -a genuine structural property of a reaction network.
At the most basic level kinetic rate equations are given by generalized mass-action functions with the partial logarithmic derivatives given by The partial logarithmic derivatives, corresponding the scaled elasticities of Metabolic Control Analysis (MCA), are restricted to constant (usually integer) values -corresponding to the kinetic order of each reaction with respect to its substrates.
We emphasize that many deterministic biochemical reaction networks can be described at the level of mass-action kinetics, such that all reaction equations comply with the functional form given above. In this case, and in the absence of mass conservation relationships, the columns of M are invariant under changes in rate constants and therefore already represent global structural properties of the reaction network.
However, sometimes the computational description of biochemical networks requires nonlinear equations -usually arising from approximations by rapid equilibrium or quasi steady state assumptions. In this case, the respective logarithmic partial derivatives are dependent on kinetic parameters and may take different values for different stationary states of the system. For example, for a generic Michaelis-Menten equation, we obtain In this case the partial logarithmic partial derivative is a non-constant quantity that explicitly depends on kinetic parameters, as well as on the stationary state of the system.
Nonetheless, our framework is also applicable in this case -as shown in the previous sections.
We note that the requirements for a given network to exhibit global concentration robustness depends to some extend on the interpretation of the logarithmic partial derivative. In particular, we can distinguish between two scenarios with respect to the interpretation of the logarithmic partial derivatives. Within the most strict assessment of global concentration robustness, we can assume that all partial derivatives are unknown and possibly variable quantities. This assumptions then also extend to generalized mass-action kinetics, such that respective logarithmic partial derivatives are not necessarily assumed to be constant quantities. In this case, the system exhibits concentration robustness also in the face of deviations from mass-action kinetics.
However, within a less stringent scenario -usually adopted within this work -we assume that the partial derivatives of mass-action rates are constant quantities that are part of the topology of the respective network. In this case, global concentration robustness may not be guaranteed for possible deviations from the assumed exponents.
We emphasize that the distinction described here does not affect or restrict the application of our framework to any actual topology -but rather it highlights that different assumptions on the nature of the rate equations may lead to different results with respect to global concentration robustness. See also Section IV for further analysis.

III. APPLICATION OF THE FORMALISM TO A SIMPLE REACTION NET-WORK
To exemplify our formalism, we consider a simple example for output invariance, as shown in Fig. 1 of the main text. Here, an output variable a is subject to slow perturbations P in its synthesis rate. Rather than fine-tuning the value of p, we look for conditions, such that an intermediate variable m compensates perturbations and ensures perfect robustness with respect to the perturbation. The pathway is described by differential equations of the form In the following, we require the system of differential equations to be well-defined, that is, all rate equations comply with basic assumptions about biochemical rate equations and the system gives rise to a positive steady state for any value of the perturbation P . Apart from these basic requirements, our method does not require to further specify the precise functional dependencies of the rate equations.
As depicted in Fig. 1 of the main text, the perturbation p acts on the rate v +a with an (unspecified) nonlinear dependency v +a = v +a (p). The elements of the perturbation vector P are defined as the logarithmic partial derivatives of the rate equations with respect to the perturbation. We use the abbreviation η := ∂ ln v +a /∂ ln p and obtain We emphasize that our analysis does not require knowledge of the precise value of η, which usually depends on the specific functional form of the rate equations, the kinetic parameters, and the strength of the perturbation.
To obtain a condition for perfect robustness of the variable a with respect to variations in p, we follow the steps described in the main text. First, we determine the basis vectors of the nullspace of N · diag(v s ). In practice any linear algebra software can be employed, such as the built-in function null(N) in matlab (The MathWorks). These basis vectors form the columns of the matrix K. For the simple example, Eq. (30), the space spanned by the columns of K is given by with Next, we consider the matrix M , with elements corresponding to the logarithmic partial derivatives of the rate equations with respect to the variable m. With the abbreviation We note that here the reactions v ±m do only depend on the output variable a and not on the variable m -a well-known prerequisite for perfect adaptation. Finally, we can assemble the invariant subspace, I = col(M ′ |K ′ ), using elementary matrix operations (EMO) A parameter-independent representation I of the invariant perturbation space I is given by In general the matrix representation I of the invariant perturbation space is not unique.
Any perturbation vector, P , that is element of the invariant perturbation space does not affect the stationary output of the system. Output invariance is tested by the rank condition, rank(P |I) = rank(I), that reads -after performing elementary matrix operations to remove the unknown functional dependencies -as The rank condition shows that the condition for global robustness of the output variable, a, is fulfilled. The example implements an integral feedback mechanism to realize a perfectly adapting reaction system. Here, the intermediate variable m acts as the 'integrator', with the crucial requirement that the rate of change of m is independent of the variable m itself. In our example, this implies that the rates v ±m do not depend on the variable m (or, equivalently, that the dependence is with the same order) -otherwise perfect robustness cannot be achieved.
Indeed, we can give a counterexample to robust behavior if we assume that the degradation In this case, Eq. (33) is modified and now reads with γ := ∂ ln v −m /∂ ln m denoting the unknown nonzero derivative. We obtain Since β and γ are both unknown and variable quantities, a parameter-independent representation I of the invariant perturbation space is restricted to two dimensions, therefore the variable a does not exhibit global robustness.

A. Notation
In this section, we provide a rigorous mathematical derivation of the conditions for local and global concentration robustness as obtained in Section II of the SI. Thereby, we consider biochemical network models with conserved moieties, and with reaction rates that may be a mixture of generalized mass action (GMA) with fixed exponents and arbitrary elements.
First, in Section IV B, we derive the previously described rank condition as sufficient and necessary condition for local concentration robustness. Then, in Section IV C, to deal with the problem of global concentration robustness, we define the invariant perturbation space such that it is independent of the network's parameters and the non-GMA part of the reaction kinetics. Based on this definition, we derive a sufficient and (in a structural sense) necessary condition for global concentration robustness.
For ease of notation, we will restrict this section to a scalar perturbation. However, this is not a restriction of generality: if multiple perturbations are present, the condition can be evaluated for each perturbation individually, and the network is robust against combined perturbations if and only if it is robust against each perturbation individually.
The image of a matrix A is denoted by im A, its kernel by ker A. A diagonal matrix with diagonal entries taken from the components of a vector a is denoted by dg a. The sum of two subspaces W 1 and W 2 of R n is defined by A biochemical network model with conserved moieties is given by the differential algebraic Thereby, x ∈ R m is the vector of independent concentrations, N ∈ R m×k the stoichiometric matrix, and v the reaction rate vector. For the conserved moieties, we have x T ∈ R n as the vector of total concentrations, while x D is the vector of dependent state variables, related to the independent state variables x via the link matrix L ∈ R n×m [3]. We consider the effect of a perturbation to the network via the variable p ∈ R, which is element of a connected and open set P ⊂ R. Note that we have explicitly accounted for the possibility that the vector of total concentrations x T may depend on the perturbation variable p. The model (42) can be transformed to an equivalent differential equation by substituting the dependent variables, The species vector x is split into output variables and intermediate variables by writing where x A ∈ R m A are the output variables and x M ∈ R m M are the intermediate variables.
According to this splitting, we introduce the matrices J A , J M ∈ R m×m given by where I m is the identity matrix of dimension m, yielding   0 We assume throughout that there exists a perturbation-dependent positive steady state in the independent variables x, given by a function x s : The steady state map for the dependent state variables is defined as For ease of notation, we definev as a shortcut for the steady state reaction rates.
Throughout this section, we assume that the steady state x s (p) is asymptotically stable for all p ∈ P. Note that this implies that the Jacobian of the right hand side of (43) in steady state, is invertible for each p ∈ P.

B. Local concentration robustness
Local concentration robustness at a perturbation value p is defined as the property that the derivative of the steady state output variables x A s with respect to the perturbation is equal to zero at p. The formal definition is given as follows.
Definition 1. The network (42) is said to have local concentration robustness at p ∈ P, if In order to derive a necessary and sufficient condition for local concentration robustness, we introduce the following matrices: With these definitions, we next introduce the invariant perturbation space, a central network characteristic for concentration robustness. In words, the invariant perturbation space is the vector space of all infinitesimal directions in which perturbations can act on the steady state reaction rates without affecting the concentration values in steady state. The formal definition makes use of the matrices defined in (52) and is given as follows.
Definition 2. The space is called the local invariant perturbation space of the network (42) at p.
We then have the following result as a condition for local concentration robustness of network (42).
Theorem 1. The network (42) has local concentration robustness at p ∈ P, if and only if Proof. First, from differentiating the steady state equation (47) with respect to p we obtain We denote Making use of the fact that J A + J M = I, we rewrite (56) as Necessity. Under the condition that J A H(p) = 0, we find that All P (p) which solve this equation are given by with a 1 ∈ ker(N dgv(p)), which implies (54).
Sufficiency. The condition (54) implies that we can write P (p) as with a 1 ∈ ker(N dgv(p)) and a 2 ∈ im J M . Thus, one particular H(p) which solves (56) is Since the Jacobian (50) is invertible, (56) has a unique solution H(p) for each P (p), which is given by (61). This proves local concentration robustness.
In the next step, we turn to the property of global concentration robustness. Essentially, a network is said to have global concentration robustness if the output variables x A in steady state are constant over the perturbation set P. In addition, we are interested in a characterization of global concentration robustness which is given by the network structure alone, and does not depend on exact parameter values and specific functions, for example Michaelis-Menten or Hill kinetics, for generic reaction rates.
To this end, we separate reaction rates into a mass action and a generic part. Thereby, the characterization of global concentration robustness should not depend on the rate constants of the mass action part nor the exact choice of mathematical function of the generic part.
However, the characterization may depend on the interaction structure of the network, i.e. which concentrations affect which reaction rate, and the stoichiometric coefficients of the mass action part entering the reaction rate expression as exponents, which are both attributed to the structure of the reaction network. Thus, in the following, the reaction rate vector v is assumed to be composed by elements of the form with i = 1, . . . , k, where Φ i represents the generic part of the reaction rate and also includes the rate constant for the mass action part, and the rest represent the concentration dependent terms in the mass action part, with stoichiometric coefficients a ij and a D ij . Before considering global concentration robustness, we first introduce the weaker notion of concentration invariance, i.e. the property that the output concentrations are constant over P, but not necessarily independent of network parameter values and choice of generic reaction rate expressions.
for all p ∈ P.
Proof. Since P is connected and x s is assumed to be a continuously differentiable function over P, the condition that x A s (p) is constant over P is equivalent to for all p ∈ P. Then, the result is a direct consequence of Theorem 1: if (63) is satisfied, Theorem 1 assures that (64) holds for any p ∈ P, thus we have global concentration invariance.
Conversely, if P (p) / ∈ I(p) for somep ∈ P, then by Theorem 1 it follows that J A x ′ s (p) = 0, and thus x A s (p) is not constant over P.
We call (63) the rank condition for global invariance, since it can be tested numerically by checking that rank(P (p)|I(p)) = rank I(p), where I(p) is a matrix whose columns span the space I(p). Concerning the mass action part of the reaction rates, the interaction structure of the network is characterized by which stoichiometric coefficients are equal to zero. For the further steps, a similar notion is also needed for the generic part of the reaction rates, and is given by the following definition.
In the following, we will frequently consider a biochemical network with the same species vectors x and x D as well as the same stoichiometric matrix N and link matrix L as the original network (42), but with potentially different reaction rates and a different vector of total concentrations. This network is described by the equatioṅ whereṽ has elements given bỹ The formal definition for global concentration robustness is as follows. In what follows, we will derive a necessary and sufficient condition for global concentration robustness based on the rank condition established above for global concentration invariance.
First, the following technical definition for structural properties of a function Φ is introduced.
Definition 6. Given a matrix ϕ ∈ R k×m and a function Φ : for all i, j.
The result on global concentration robustness uses the two matrices M and P defined as follows, which depend on matrices ϕ x ∈ R k×m , ϕ x D ∈ R k×n , ϕ x T ∈ R n , ϕ p ∈ R k , and vectors Proof. Sufficiency: Consider a network given by the equations (66) with reaction rates as in (67), whereΦ andx T are structurally equivalent to Φ and x T , respectively. Let x s (p) and x D s (p) be steady state concentrations of this network andv(p) the corresponding steady state reaction rates. DefineΦ(p) =Φ(x s (p), x D s (p), p). Also define the following matrices: and note that ϕ p ∈ S( ∂Φ ∂p ), ϕ x ∈ S( ∂Φ ∂x ), ϕ x D ∈ S( ∂Φ ∂x D ), and ϕ x T ∈ S( ∂x T ∂p ). The matrices P , Q, D, and M from (52) for the network (66) are computed as follows: The condition P (ϕ x D , ϕ x T , ϕ p ) ∈ im M (ϕ x , ϕ x D , x D , x) + ker(N dg α) implies that P (p) ∈ im M (p) + ker(N dgv(p)), for all p ∈ P, implying that the network (66) has global concentration invariance. Thus, by Definition 5, the network (42) has global concentration robustness.
Necessity: Assume that the condition P (ϕ , and ϕ p ∈ S( ∂Φ ∂p ): Next, consider the network (66) with reaction rates as in (67), whereΦ is chosen as Choose anyp ∈ P and let for i = 1, . . . , k. Furthermore, definē and letx for i = 1, . . . , n. Note thatΦ and Φ as well asx T and x T are structurally equivalent due to the structural constraints on the matrices ϕ x , ϕ x T , ϕ p , and ϕ x T . With the definition ofk i in (76),ṽ(x,x D ,p) = α. Then, since α ∈ ker N , and alsox T (p) = Lx +x D , the network (66) withΦ,x T as just defined has a steady state given by and steady state reaction ratesv It remains to show that (66) does not have local concentration robustness at the perturbationp. The matrices P and M defined in (52) are computed for the network (66) as follows: Thus, from condition (74), we find that Note that the condition (70) can again be rephrased as the rank condition rank(P |I) = rank I, where I is a matrix whose columns span the space im M + ker(N dg α).
From the stationary equations of the reaction rates, given in Materials and Methods of the main text, we determine the logarithmic expansion coefficients for the state variables M and the perturbations P . In the following only nonzero coefficients are considered and the results are summarized in Table I.
We make use of Kronecker's delta, defined by δ ij = 1 for i = j and zero otherwise. If perturbations act on individual components we arrive at which results in the matrix shown in Fig. 2 of the main text.
The non-zero entries of M are given by  (2) ). Note that the rate equations contain nonlinear functions k A = k A (ATP), P s = P s (m, L) and A c = A c (A T , W T , T ). Lowercase Greek letters denote the logarithmic expansion coefficients that arise from (unspecified) nonlinear dependencies. We emphasize that, though in this case the precise functional dependencies are known, our framework does not require to specify the exact functional form of the dependencies.
As demonstrated in the main text, the specific topological organization, with kinetic dependencies summarized in Table I, allows the output of the chemotactic pathway to be robust against diverse perturbations that would otherwise impede the functionality of the network.
As can easily be ascertained in Table I each perturbation corresponding to the columns of P (1) is an element of the invariant perturbation space. In contrast to this, perturbations corresponding to the columns of P (2) are not within the invariant perturbation space, hence the pathway is not robust against fluctuations in these components. However, the organization of the pathway ensures that the pathway is indeed robust against concerted fluctuations in these components. See main text for details.
We emphasize that the uncovered design principle of the chemotaxis pathway is in contrast to the more straightforward possibility to utilize an extensive cellular machinery to 'finetune' quantities that appear as parameters in the equations, such as protein concentrations or ATP availability. We further note that the robustness requires some logarithmic derivatives to attain specific values, as encoded for example by mass-action kinetics, while other functional dependencies may remain unspecified.

B. Growth conditions
All strains were grown under standard chemotaxis conditions [11,13] at 34C in a rotary shaker to mid-exponential phase (OD600 ≈ 0.48) in tryptone broth (TB) supplemented with 100µg/ml ampicillin and indicated amounts of IPTG.

C. FRET measurements
Cell preparation, FRET measurements and evaluation of FRET data were performed as described previously [12,14] on a custom-modified Zeiss Axiovert 200 microscope.

D. Quantification of gene expression
Expression of fluorescent reporter proteins in individual cells was quantified as described before [5] using fluorescence imaging on an AxioImager fluorescence microscope equipped with an ORCA AG CCD camera (Hamamatsu).   [5]. The resulting rescaled kinase activity is shown in Fig. 4B.

A. Two-component systems and implications for synthetic biology
One of the merits of our approach is to guide the design of perfectly robust signaling circuits -with important implications for synthetic biology. To exemplify the construction of robust signalling networks, we briefly consider instances of two-component signal transduction systems. Bacterial two-component systems typically consist of a membrane-bound sensor kinase that senses a specific stimulus and a cognate response regulator that modulates the signal response. The robustness of individual bacterial two-component system with respect to concentration fluctuations was investigated previously [2,9].
Recently, Skerker et al. [10] described a method that allows for the rational rewiring of the specificity of two-component systems. Such rational design of the output responses of two-component systems is a major step forward in the design of protein-based synthetic pathways, with exciting potential applications in synthetic biology and biotechnology [4].
However, the rational design of two-component systems will also necessitate to engineer robustness of the rewired pathways with respect to possible detrimental fluctuations -taking into account that synthetic circuits are not a product of evolution. In this respect, of particular importance are fluctuations in compound concentrations that arise from stochastic variations in transcription and translation, as well as from other sources, such as variations in division. In fact, it seems highly desirable to implement any altered topology such that the expression of the individual proteins has no effect on the (rationally designed) inputoutput relationship of the network. In this case, the only requirement for the method of Skerker et al. [10] to generate perfectly robustness networks is a sufficiently high expression of any of the involved proteins -without the need to fine-tune any of the precise expression levels.
Our framework is able to straightforwardly account for the mechanisms of robustness of such (networks of) two-component systems. In Fig. 1 two variants  (but see below for the full solution), the system is described by the two differential equations for the independent state variables, H P and R P , and the mass conservation relationship Both topologies only differ in kinetic dependencies of the rate equations, specifically for the systems shown in Figs. 1A and 1B, respectively. The right nullspace of the scaled stoichiometry is identical for both systems and can be represented by a matrix K that solely consist of the 1-vector (here the column has already been normalized to remove dependencies on the stationary flux distribution).
We aim to test for invariance of the pathway output -the phosphorylated response regulator (R P ) -with respect to variations in total component concentrations H T and R T . We proceed as described in Sections II.B and II.C.
We start with the topology shown in Fig. 1A. In this case, the matrices of logarithmic partial derivatives are given as resulting in a matrix M defined as (see Section II.B) where α = −H p /H. The perturbation vectors with respect to the total concentrations H T and R T are given as with γ H := ∂ ln ν 1 /∂ ln H T and γ R : The rank condition for global concentration robustness with respect to the total concentration of the histidine kinase (H T ) therefore reads Obviously, the equation cannot be fulfilled for arbitrary values of γ H and α, hence the system shown in Fig. 1A does not exhibit global concentration robustness with respect to the total concentration of the histidine kinase (H T ).
The same conclusion can be reached for the total concentration of the response regulator (R T ). The rank condition for global concentration robustness reads Again, the equation cannot be fulfilled for arbitrary values of γ R and α, hence global concentration robustness is not achieved.
A different scenario is shown in Fig. 1B. Here, a bifunctional histidine kinase implies that the unphosphorylated kinase acts as a phosphorylase for the response regulator. Repeating the calculations shown above, we obtain resulting in, and with definitions given above. We emphasize that indeed γ H : Testing the rank condition for both perturbations, reveals that thus both conditions are fulfilled for any non-zero values of γ H , γ R , and α. The system shown in Fig. 1B indeed exhibits perfect concentration robustness of the output variable R p with respect to variations in both total concentrations H T and R T . We note that within this example the bifunctionality of the histidine kinase is crucial to achieve robustness of the pathway output -a mechanism that is functionally similar to the concerted expression of proteins adjacent on an operon observed for the E. coli chemotaxis pathway.
Elaborating on the simple system discussed above, our framework is also straightforwardly applicable to the full system, including explicit complex formation. We consider the three processes H → H P (109) corresponding to system of 6 variables and 7 reaction rates. The system of differential equations is given as and supplemented by the mass conservation relationship Here the concentrations of unphosphorylated components, R and H, are chosen as dependent variables. We note that the choice of dependent variables is not unique. Alternative choices lead to identical results, provided that the dependent variables are chosen such that the matrix M D is of maximal possible rank.
To test for robustness of the pathway output, we first evaluate the nullspace of the stoichiometry, where the representation has already been normalized and dependencies on the stationary flux distribution were already removed. The logarithmic expansion coefficients for the state variables M D are summarized in Table II. Instead of computing the expression for M = M D − D · L ′′ , we utilize a direct approach to judge output robustness of the network. Since we are only interested in the vector space spanned by the matrices M and K ′ , and not in a particular representation, we note that the expression for M is not required under the condition rank(D|M D |K ′ ) = rank(M D |K ′ ). As can be verified in Table II this  This result can be generalized to a generic strategy towards perfect output robustness for engineered protein networks. Utilizing rewiring of substrate specificity [10], in addition  TABLE II: The logarithmic expansion coefficients with respect to the dependent and independent state variables, D and M D , respectively, along with the nullspace K ′ .
to ensure global concentration robustness against all total concentrations within a signaling network. This simple condition allows to establish whether a rewired network will exhibit the desired functionality without the need to fine-tune expression levels. Our framework is able to guide the necessary network extensions to guarantee robust network functionality.

B. Conservation relationships and robustness of mass-action systems
A particular application of our framework relates to global concentration robustness of mass-action systems with respect to total conserved concentrations, as recently also discussed elsewhere [8]. Our framework allows to derive a simple principle that allows to judge for global concentration robustness and is able to guide the necessary network extensions to design perfectly robust networks. We consider a system as described in Eq. (7) in Sec- and aim to test for robustness of the output variables x A with respect to perturbation in the total concentration of each molecular component x T = (x T 1 , x T 2 , ..., x T n ). For simplicity, we assume that all rate equations are given by generalized mass-action (GMA) kinetics, Consequently, the partial logarithmic derivatives with respect to dependent and independent where K ′ denotes a largest parameter-independent representation of the nullspace, the system exhibits global concentration robustness with respect to perturbations in x T . The reason is that any perturbation in total concentrations can be represented by a perturbation in the dependent variables x D . Since D is already an element of the invariant perturbation space, due to the condition Eq. (119), any such perturbation is necessarily also an element of the invariant space, hence global concentration robustness with respect to total conserved concentrations is guaranteed.
As compared to alternative methods [8], our approach has the advantages that it is (i) conceptually considerably simpler, (ii) numerically straightforward to test by standard methods of linear algebra (rank conditions), and (iii) straightforwardly guiding modification of the matrix D to ensure perfect concentration robustness.

C. Complex perturbations and temperature compensation
One of the merits of our approach is that it is not restricted to a particular type of perturbation, but is applicable to a wide variety of detrimental influence that potentially impede network functionality. While our focus is mainly on variations in native expression levels -as one of the dominant sources of variability in living cells -our framework also accounts for any other perturbation that can be expressed in terms of the logarithmic expansion coefficient of the rate equations.
Relevant applications include the retroactivity of signaling circuits [15] as well as detrimental pathway crosstalk [1,6]. Both issues also relate to scenarios where signaling pathways utilize common resources, such as ATP to provide energy. Within our framework, any such possibly detrimental influence can be considered as a perturbation -allowing the identification or construction of an appropriate topology that compensates for the corresponding perturbation.
A particular intriguing example of a complex perturbation is given by variations in temperature. A change in temperature usually affects all rate constants simultaneously -making the prediction of perfectly robust topologies a difficult task [7]. In the simplest case, we may assume that each reaction rate follows the Arrhenius equation, that is, the temperature dependence of a reaction rate v i can be described by a multiplicative factor where E i denotes the activation energy for the ith reaction, A i a proportionality constant, R the gas constant, and T the absolute temperature in Kelvin). Here, the activation energy is a constant for each reaction that does not further dependent on the temperature or the stationary state. In this case, we can straightforwardly construct the perturbation vector P T with respect to changes in temperature, with elements Specifically, each element of P T explicitly depends on the temperature. However, the direction of the vector P T , with for a system consisting of r reaction rates, does only depend on the (constant) activation energies of the reaction. Hence, using our condition for output robustness it is straightforward to judge global temperature compensation of a biochemical network.
We note that in practise a straightforward application of the Arrhenius equation is often not appropriate. In this case, for arbitrary dependencies k i (T ) our framework is still applicablejust as for any other complex perturbation that acts on many reaction rates simultaneously.

D. Robustness of bi-and multistable systems
As yet, the focus of our approach has been networks that exhibit a globally stable stationary state, characterized by x s and ν s , for all parameters. However, many signaling networks exhibit bi -or multistable dynamics, giving rise to two or more stationary states.
As a particular merit, our approach is still applicable in these situations -without requiring substantial modifications. In particular, global uniqueness of the stationary state was not a necessary precondition to derive the requirement for global robustness, hence the condition for global concentration robustness applies to any locally stable state that fulfills the steady state condition.
However, our framework does rely on invertibility of the Jacobian matrix to ensure that a gradual change of intermediate variables may not affect the set of designated output variables. This does not hold in a situation in which the system, under the action of a perturbation, undergoes a bifurcation that results in a non-invertible Jacobian (as, for example, a saddle-node bifurcation). With respect to the application on multistable systems, we therefore have to introduce the additional constraint that all perturbations must be such that the transient response in systems variables after the perturbation remains within the basin of attraction of the respective state. Once the perturbation is sufficiently strong to allow the system to cross the attractor boundary, the robustness of the state is lost. However, in this case the system usually adopts another stationary state -which again exhibits perfect concentration robustness with respect to perturbations of large magnitude. In this sense, we are in the favorable situation that our framework allows to construct robust bi-or multistable systems that are still capable to switch between states. Also, since we are mainly concerned with slow perturbations with respect to the intrinsic timescales of the system, the robustness of each stationary state is maintained even for perturbations of comparatively large magnitude -provided the variables remain within the basin of attraction of the respective state.
As a guideline for the construction of robust bi-or multistable systems, we further note that such systems usually involve strong nonlinearities. Here, it is of considerable advantage to utilize a robust ("output") variable within the feedback mechanism. In this way, the requirement to fine-tune a highly nonlinear feedback is circumvented.
We illustrate the design of a multistable robust system using a simple example based on a generic two-component system. We again consider the system of differential equations discussed above The vector of kinetic rate equations ν is analogous to the robust topology shown in Fig. 1B.
As the only modification, we assume that the active response regulator R p is capable to sensitize the receptor, thereby enhancing the phosphorylation of H. This modification does not impede the robustness of the system with respect to changes in the total conserved protein concentrations. In particular, the derivative towards the variable R p does not enter the invariant perturbation space or its construction, hence the invariant perturbation space is identical to the scenario discussed above.
Solving for the steady state of the system using the relationship ν 1 = ν 3 at steady state, we obtain With a particular choice of f (R p ), for example, and n = 4 we obtain and the system indeed exhibits bistable dynamics that is independent of the total concentrations H T and R T . A bifurcation diagram is shown in Fig. 2.
The system allows for several conclusions: (i) Concentration robustness holds for any stationary state of the system, including unstable states. (ii) Strong nonlinearities should be confined to robust variables. As these variables do not enter the construction of the invariant perturbation space, the need for fine-tuning the respective parameters is circumvented. (iii) Robustness of the system is lost at the bifurcation. In this case, the state looses stability and adopts another state, which again exhibits non-local robustness. Hence, robustness is restricted to perturbations within the respective basin of attraction. (iv) While these conclusions may be obvious for the simple example discussed above, our reasoning likewise applies to systems of large size and is amendable to an algorithmic solution. For many cases this breakdown can be made more explicit, when accounting for conservation equations. For example the active form of a protein y 1 is related to the total concentration by y T = y + y 1 , with y the inactive form. Thus, we have y 1 ≤ y T . As the rank condition only uses differential forms, the latter inequality has to be additionally fulfilled.
It is further important to note that robustness as defined by the rank condition

B. Fine-tuned global robustness
For some reaction networks it might not be clear beforehand if the parameter dependent logarithmic derivatives -denoted by Greek symbols -are indeed independent. In this case there could exist a functional dependence between some of the non-constant elements in the matrix (P |M |K) that increases the invariant subspace. An example is given by the following (rather artificial) system that includes a specific perturbation parameter, 0 < p < 1, In the stationary state we can substitute x s = − ln(p) into the second equation get y s = k 1 /k 2 . Thus the systems output is obviously independent of p. Although this systems shows global robustness against p, the reaction network involves fine tuning, in the sense that two independent reactions share the same rate constant, k * . Using our formalism in a straight forward manner, we arrive, after introducing the independent functions α and β, at However, as α := x s , β := [ln p] −1 and x s = − ln p, we see that α = β −1 . After multiplication of the P -column with β −1 we see that the rank condition is indeed fulfilled. This fine tuned global robustness can be numerically detected by inserting the actual values for the Greek symbols for one stationary state of the system. If robustness with respect to P differs in this case to randomly chosen Greek symbols, then there exists a hidden dependency among the parameter dependent logarithmic derivatives. This example shows that global robustness of reaction network can also emerge from a combination of fine-tuning and network structure.
Here the hyperplane spanned by the column vectors of (M |K) is rotating in space and the perturbation, P , is such that it constantly follows the rotation of this hyperplane for all p.
In general, our assumption of independence between the partial derivatives represents a worst-case scenario and is sufficient for global robustness.

C. Numerical test of the rank condition
In practice, a simple numerical test for the rank condition Eq. (6) can be performed by treating the independently varying logarithmic derivatives as random variables. In the case of the example in Section II B this can be realized by redefining the Greek symbols as random variables that take values according to a uniform distribution U (−1, 1). In doing so, it has to be guaranteed that each randomly chosen value must be sufficiently different from zero and sufficiently different from all other random values. The minimum differences between the random values are set by the numerical precision and avoid that by chance almost equal values are assigned to linear independent logarithmic derivatives.
D. Determination of the parameter-independent nullspace K ′ So far, the nullspace K -determined by N · diag(v s ) · K = 0 -depends on the particular flux distribution v s . To identify the parameter-free conditions for structural robustness, it is desirable to identify the subspace of K that is independent of the particular state v s .
To this end, we first construct the nullspace K from the nullspace of the original stoichiometric matrix, N · K N = 0, using the transformation where q denotes an invertible (k − m) × (k − m) matrix, corresponding to a basis transformation of the non-unique representation of the nullspace.
At this point it may be argued that the flux dependency of K can be partially removed by using elementary matrix operations. However, it is far to restrictive to treat the stationary fluxes v s = (v s 1 , ..., v s k ) as unknowns as the fluxes are not linear independent v s = K N · α , with α the elementary flux coefficients and dim α = dim v s − rank(N ) = k − m. We therefore utilize the freedom of constructing the matrix q in order to generate the largest possible subspace K (1) that is independent of reaction fluxes. We first rewrite Eq. (131) as The parameter independent subspace of K defines K (1) and requires that the elements of K (1) are independent of α. This implies the condition ∂ αs K (1) nm = 0 = ∂ αs j K N nj q jm j K N nj α j = j K N nj (∂ αs q jm ) for all indices s, which can be rewritten to give j K N nj ∂ αs q jm = j K N nj q jm j K N nj α j const K N ns (135) that in turn takes -by defining K N j as the j-th column of K N -the alternative form The latter statement requires q jm ∝ α j for all m and thus at least one column vector of q exists, given by q * = α, that results in a parameter free representation of K. By Eq. (133) we obtain the first parameter free nullspace vector, K (1) n1 = j K N nj α j / j K N nj α j = 1 for all rows, n. A k-dimensional vector consisting only of ones is thus always part of the invariant subspace and reflects an obvious invariance property of all stationary networks: multiplication of all fluxes in the network by the same factor does not change any stationary state variable of the network.
As this invariance property can hold also locally, we can separate q * in the maximum number of orthogonal column vectors q * 1 = (α 1 , ..., α l 1 , 0, ..., 0) T , q * 2 = (0, ..., 0, α (l 1 +1) , ..., α l 2 , 0, ..., 0) T , ..., q * z = (0, ..., 0, α l z−1 +1 , ..., α k−m ) T such that Eq.(136) holds. The independent columns {q * 1 , ..., q * z } can be determined by a block matrix representation of K N which is obtained by resorting columns and rows such that the resulting matrix has only zero entries to all sides of each block. Note that the blocks are in general not square. This resorting leads to Thus the column vectors of the matrix K (1) that indicate robustness to a fold change in several fluxes can be constructed by the surprisingly simple transformation where the 1's denote a vector of ones with dimension set to number of rows of the corresponding block. K (1) . Next we complete K (1) and identify the complementary nullspace, K (2) , to K (1) such that both spaces together span the complete nullspace of N · diag(v s ), K = (K (1) |K (2) ).
We illustrate the construction of K (1) and K (2) by an example of a reaction network which gives rise to two alternative stationary flux distributions. The example consists of two The stoichiometric matrix and its corresponding nullspace are given by showing explicitly the two above mentioned flux distributions. The nullspace K (1) consists of one column vector K (1) = (1 1 1 1 1 1) T as obviously K N can be grouped only in one single block. As K = [diag(v s )] −1 · K N the i-th row of K N is weighted by (v s i ) −1 = (K N · α) −1 i . As the α i can vary independently under perturbations, the identity v s i = v s j holds if the i-th and j-th row of K N are identical. Identical rows of K N thus weight the the rows of K the same inverse flux, whereas all other rows of K carry different weights. As the dimension of the joint vector space spanned by K (1) and K (2) must have the same dimension as Kwhich by construction has the the dimension of K N -we can choose only one column of K N to construct K (2) in this example. Taking the second column and indicating the unknown inverse fluxes by Greek symbols we obtain with β 1 = (v s 1 ) −1 = (v s 4 ) −1 , and β 2 = (v s 5 ) −1 = (v s 6 ) −1 . If we multiply the first column of K by β 2 and subtract the second column we obtain β 2 K 1 − K 2 = (β ′ 1 β 2 β 2 β ′ 1 0 0) T with β ′ 1 = β 2 − β 1 . This result would have been obtained by taking the first column of K N to construct K (2) . We note that K (2) is in general not part of the invariant subspace, due to its dependence on the (local) flux distribution. This dependence is obvious in the example above, where the two column vectors of K span a two dimensional hyperplane and one vector cause the hyperplane to rotate in six dimensional space under a change in stationary flux distributions. The method introduced so far allows to construct in a systematic way the nullspace K such that the dependence on stationary fluxes -that is the number of unknowns β i -is reduced to a minimum. The subspace spanned by the columns K can be further reduced by elementary matrix operations (EMO) with δ = β 2 /β 1 .