A Network Characteristic That Correlates Environmental and Genetic Robustness

As scientific advances in perturbing biological systems and technological advances in data acquisition allow the large-scale quantitative analysis of biological function, the robustness of organisms to both transient environmental stresses and inter-generational genetic changes is a fundamental impediment to the identifiability of mathematical models of these functions. An approach to overcoming this impediment is to reduce the space of possible models to take into account both types of robustness. However, the relationship between the two is still controversial. This work uncovers a network characteristic, transient responsiveness, for a specific function that correlates environmental imperturbability and genetic robustness. We test this characteristic extensively for dynamic networks of ordinary differential equations ranging up to 30 interacting nodes and find that there is a power-law relating environmental imperturbability and genetic robustness that tends to linearity as the number of nodes increases. Using our methods, we refine the classification of known 3-node motifs in terms of their environmental and genetic robustness. We demonstrate our approach by applying it to the chemotaxis signaling network. In particular, we investigate plausible models for the role of CheV protein in biochemical adaptation via a phosphorylation pathway, testing modifications that could improve the robustness of the system to environmental and/or genetic perturbation.


Introduction
Biological systems in general show various types and degrees of robustness to environmental changes, meaning that they continue to function even when changes in the environment occur. This imperturbability is often accompanied by robustness to genetic perturbations, meaning that progeny function even though their genotype is not identical to the parent genotype [1][2][3][4]. Both features play an important role in evolutionary biology. While the former is a direct outcome of selection, the relationship between evolution and genetic robustness is likely to be indirect for low functional mutation rates [5][6][7] since selection acts only on the phenotype of an organism and not its genotype [8].
In this study, we develop a computational experiment to investigate the plausibility of this hypothesis, that there is a general correlation between environmental and genetic robustness, and provide a quantitative measure of the degree of correlation, if any. In more detail, we shall show that the presence of a specific dynamic network characteristic in networks is associated with a better correlation between genetic and environmental robustness than found in networks where it is absent. Rather than focusing on a particular system in a specific organism, we choose one function of interest: The ability to attain steady state output for constant input. If a network capable of carrying out this function is robust to external environmental perturbations, what is the probability that it is also robust to internal (e.g., genetic) disruption? To be specific, we define environmental robustness of a biological network as the ability to maintain an output in the face of input perturbations. Genetic robustness is defined as the ability of a biochemical system to maintain the same output in the face of genetic mutations represented as rate constant changes in the equations representing it. This representation of a mutation as a jump from one set of parameters to another is a standard assumption [29].
For mathematical convenience, we restrict our discussion to Michaelis-Menten type networks as they are likely to reach a steady state under constant inputs relative to general networks without sigmoidal saturation. Such networks were also used in the analysis of three node biochemically adaptable networks by Ma et al. [30]. The sensitivity of biochemical kinetic models to parameter perturbations has been intensively investigated [29,[31][32][33][34][35]] as a mathematical model of a biological system should be able to reproduce the function of interest or fit experimental data with a minimal need for parameter fine-tuning [35,36]. Systems of biochemical adaptation [30,[37][38][39][40] have been of interest in particular.
Defining a topology to be a graph of interactions independent of parameter values, we test a large number of random N-node topologies for networks capable of reaching a steady state both under constant input concentrations and after a persistent step change in these input concentrations. We define a network as a topology with a specific set of parameters. Each network is given a numerical value for its level of robustness to input and parameter perturbations. The level of robustness of the topology is determined by averaging over this value obtained from its corresponding networks. In particular, we differentiate between networks that show a transient response to a step change in input and those that do not. We find that there is a statistically significant model II regression between the level of robustness to input of a topology and its level of robustness to parameter perturbations that has a steeper slope in networks with a transient response. Our results may be relevant to the discussion about the relationship between the need to survive in a constantly changing environment and the evolution of genetic robustness.
There is a large literature on functional motifs that are necessary for a biological system to carry out specific tasks [30,[41][42][43][44][45][46][47][48][49]. Here, we test all possible 3-node topologies to find the particular motifs that are of use in achieving both robustness to input and parameters. Having established the correlation between environmental and genetic robustness, we ask if there are topologies sharing certain sets of motifs/architectures that show stronger correlations than others. Ma et al [30] computationally explored all possible topologies of 3-node Michaelis-Menten enzymatic networks for motifs that can best accomplish biochemical adaptation. Using our results on this correlation between different sets of architectures we refine the list of motifs of biochemical adaptations previously published [30].
Our approach can be used to select/reject plausible/improbable models of a system of interest. We demonstrate this via a comparative study of bacterial chemotaxis signaling systems. Chemotaxis is a process generally used by bacteria to sense changes in their chemical environment [4,[17][18][19][20][21][22][23][24]. Chemotactic signaling is a well-studied system, but most of the focus has been on the chemotaxis network of the Escherichia coli (E. coli) bacterium [4,[17][18][19][20] despite the fact that chemotactic signaling pathways differ between species [21][22][23][24]. For instance, CheV is a chemotaxis protein found in many bacteria but not in E. coli. In many species, it was shown that CheV, or a variant of it, plays a role in biochemical adaptation during chemotaxis via its phosphorylatable receiver domain [24,50,51]. However, the exact mechanism is still not known [24]. Here, we compare the coarsegrained network of E. coli chemotaxis with several others involving CheV phosphorylation. We draw conclusions based on the resultant values of robustness to both input and parameter perturbations and the correlation between them.
In summary, we provide extensive evidence for a mathematical principle stating that, statistically speaking, dynamical systems that are biochemically adaptable are also genetically robust. We apply this knowledge to search for topological categories and subcategories within 3-node networks that show a particularly strong correlation and a linear relationship between their robustness to input and to parameter perturbations, and to shed more light on the chemotactic signaling pathways in bacteria. This method of searching for motifs can be extended to other functions and to bigger networks in order to find motifs that combine more complex functions necessitating larger numbers of nodes.

Results
In the current work, we sample over 50,000 topologies each of 5-node, 10-node, 15-node, and 30-node networks, and over all 3 9 possible topologies of 3-node networks. For each topology j, we average over a large number of randomly chosen parameter sets. The parameters are chosen from a uniform distribution within fixed ranges as described in the Methods section. For each network defined by topology j and parameter set P I , we compute two values: E I j,P I which is a measure of the robustness of the network to a persistent step change in input, and E p j,P I which is a measure of the robustness of the network to perturbations in its set of parameters P I . We take the geometric averages of E I j,P I and E p j,P I over the whole parameter space as a quantitative evaluation of the robustness of topology j to a step change in input and to parameter perturbations respectively ( Fig. 1).
In previous work on biochemical adaptability, it was assumed that networks that quickly respond to input change are better adapted than those with slower response [30,37]. In this work, we take a qualitative approach to avoid a bias towards larger or faster transients. A biochemically adaptable network is defined as one that is both robust to input perturbations and has a transient response to a persistent step change in input, independent of the magnitude of the transient. A step change in input ( Fig. 2A) induces three possible responses from the output dynamics (assuming a steady state can be reached): No response (Fig. 2B), a monotonic response (Fig. 2C), or a transient response (Fig. 2D). Due to possible computational noise in the time-course of the output concentration, we need an objective way to distinguish between a network with a small transient and a non-responsive or monotonically responsive one (e.g., Fig. 2E,F). To this end, we evaluate the Pearson shape correlation between the network's time-course and two model time-courses representing the dynamics of a network characterized by perfect biochemical adaptation (the red time-courses in Fig. 2E-H) and that of a monotonically responsive one (the green time-courses in Fig. 2E-H). Here, the

Author Summary
Advances in the ways that living systems can be perturbed in order to study how they function and sharp reductions in the cost of computer resources have allowed the collection of large amounts of data. The aim of biological system modeling is to analyze this data in order to pin down the precise interactions of molecules that underlie the observed functions. This is made difficult due to two features of biological systems: (1) Living things do not show an appreciable loss of function across large ranges of environmental factors. (2) Their function is inherited from parent to child more or less unchanged in spite of random mutations in genetic sequences. We find that these two features are more correlated in a specific subset of networks and show how to use this observation to find networks in which these two features appear together. Working within this smaller space of networks may make it easier to find suitable underlying models from data.
time-course represents the dynamics of the concentration of the output node from one steady state (before the change in input concentration) to a new steady state (after a persistent change in input concentration). In summary, we use the term ''transientlyresponsive'' (TR) for a network that responds to a persistent step change in input and then returns to a new steady state different than the peak response regardless of whether it is also robust (Fig. 2H) or not (Fig. 2G). Any network that does not pass the Pearson test, whether it shows no response or a monotonic one, is termed Non-Pearson (NP). A perfectly biochemically adaptable network is one that is both transiently-responsive and perfectly robust to input perturbations (Fig. 2H).

Correlation between Robustness to Input and Parameter Perturbations
We define and derive (see Methods) two quantitative measures of input and parameter robustness for each topology j: G TR For topologies with more than 3 nodes we sample over at least 50000 different ones of each size (5, 10, 15, and 30-node topologies) while 3-node topologies are exhaustively sampled. The different topologies (that have more than 3-nodes) are sampled randomly as described in Selection Criteria in the Methods section. We reject topologies with a low fraction of TR networks (F succ v2:3%, where F succ is the ratio of the number of TR networks to the total number of networks) and exclude them from any further analysis. We chose 2.3% to be the cutoff on F succ as it is the minimal value of F succ that removes clusters and outliers (Fig.  S1). With this, we are left with 2445, 7847, 18300, 19264, and 16589 3-node, 5-node, 10-node, 15-node, and 30-node topologies respectively. The networks sampled from each of these topologies are qualified as TR (Fig. 3) or NP (Fig. 4) and separated accordingly.
We find that over the parameter space of a topology j, E I j,P I and E p j,P I can span a wide range of values. Within both TR and NP networks, we find a significant (p>0.0) linear correlation between log 10 G I j ð Þ ð Þ and log 10 G p j ð Þ À Á . A comparison of the slope of the linear regression (using model II regression, in particular the ordinary least square bisector method described in [52]) shows a clear and systematic pattern between topologies of different sizes and TR and NP networks of the same size. We find that the slope increases as the size of the network increases: slope = 0.6260.09 for 3-node (Fig. 3A), 0.6460.05 for 5-node (Fig. 3B), 0.7660.04 for 10-node (Fig. 3C) 0.8060.05 for 15-node (Fig. 3D), and 0.9660.13 for 30-node topologies (Fig. 3E). The marginal error is taken as the 95% confidence interval where the variance of the slope is calculated using its estimate for OLS bisector regression derived by Isobe et al [52]. Similarly, for NP networks we obtain: slope = 0.5060.04 for 3-node (Fig. 4A), 0.4460.02 for 5-node (Fig. 4B), 0.5460.02 for 10-node (Fig. 4C), 0.5960.02 for 15-node (Fig. 4D), and 0.7060.03 for 30-node topologies (Fig. 4E). As above, the marginal error here is taken as the 95% confidence interval. The confidence intervals show that, for all sizes, the values of slopes for TR networks are consistently higher than that for NP networks of the same size and that the difference between the two slopes is significant.
The values of the Pearson correlations within TR and NP networks show no clear pattern. This is mainly due to the variability introduced by parameters whose robustness stays invariant and reducible topologies within N.3 N-node networks (see Text S1 and Discussion). Due to these caveats we are cautious about drawing conclusions based on the values of the Pearson correlation.

3-Node Correlations within TR Topologies and Corresponding Motifs
Sampling over all 3 9 possible topologies, our results show only 4153 topologies have associated TR networks. Within these topologies we find a significant linear correlation before (Fig. S2 A) and after (Fig. 3 A) introducing the cutoff, as discussed in the previous section. In what follows we show how we can extract motifs (i.e., basic topologies that may be more likely to appear in biological systems) by examining the slope of the linear regression between log 10 G TR I À Á and log 10 G TR p . Here, we show that motifs can be extracted from topologies representing the basic backbones shared by a set of topologies showing a stronger relation between environmental and genetic robustness as follows.
We first consider two known motifs, the incoherent feedforward motif (IFF) and the negative feedback loop motif (NFL) and examine their corresponding relations. IFF (Fig. 5A) is a topology wherein the output node is affected by the input receiving node via two paths, one direct and the other indirect, such that, collectively, one path is activating and the other is deactivating. This implies four subcategories denoted IFF1-IFF4. NFL (Fig. 5A) is a topology wherein a node i is activated/deactivated by another node j, and node j is deactivated/activated back by node i either directly (NFL1, NFL2) or indirectly (NFL3-NFL10). We find that the majority of TR topologies have IFF, NFL, or both IFF and NFL motifs. Only a few have neither IFF nor NFL; these are, however, robust to neither input nor parameter perturbations (Fig. 6A) and they all have low fractions of successful trials, indicating that TR networks generated from these topologies are sparse. All 4 subcategories of IFF are fairly robust to both input and parameter perturbations (results not shown). Though NFL only topologies are generally less robust than those containing IFF, a small group of them (green cluster at the bottom left of Fig. 6A) have low numbers of successful trials but are highly robust within their small TR space. When not coexisting with other robust motifs, only 4 out of the 10 categories of NFL (NFL1, NFL2, NFL4, and NFL6) are robust to both input and parameter perturbations (Fig. 6B). Seeing how NFL1 topologies show separate groups in Fig. 6B, we examine the distribution of all topologies containing NFL1 according to its 8 types (Fig. 5B). We find that NFL1 type1 topologies are strongly correlated (r = 0.92, p>0) while NFL1 type2 show separate clustering Figure 1. Flowchart for our methodology. N t = the total number of tested random topologies.Ñ N s = the total number of trials, i.e., different sets of parameters tested. * (1) for the 3-node networks the topologies are not randomly generated, rather sequentially in order to test all 3 9 possible combinations. * (2) if the network takes too long to reach equilibrium or the Jacobean matrix J I I Figure 2. Illustrating the criterion for selecting for transiently responsive networks. The Pearson test is performed to decide whether a network (a particular choice of a set of parameters for a particular topology) shows a transient response to a step change in input or not. The starting point is at steady state under a constant input concentration I 1 . At time t c the input concentration is changed from I 1 to I 2 via a step function (A). Consequently the output will either (B) not sense the change and maintain the same steady state O 2~O1 , (C) change monotonically to a new steady state O 2 =O 1 , or (D) show a transient response followed by a relaxation to a new steady state O 2 that might or might not be equal to the pre-step change steady state O 1 . A network passes the test only if it is transiently responsive. (E-H) F ad (red) is a perfectly biochemically adaptable function and F nad (green) is a monotonically changing function. If The Pearson shape correlation between the computed time course (blue) and F ad , r 1 , is bigger than that between it and F nad , r 2 , then the test is passed (E, G, H) and the network is termed a Transiently-Responsive (TR) network, otherwise the test is failed (F) and the network is termed Non-Pearson (NP). In (E) and (F), we show two cases where a biased visual inspection would deem both networks looking similar, but in reality they are distinguished by the Pearson test. A biochemically adaptable network is one that is both transiently responsive and returns to a new steady state very close to the pre-step change steady state. For example, (G) is not biochemically adaptable as there is a large difference between the pre-step change steady state and the new one. On the other hand, (H) is perfectly biochemically adaptable as it is both transiently responsive and perfectly robust, i.e. it returns exactly to its pre-step change steady state. doi:10.1371/journal.pcbi.1003474.g002 Figure 3. Correlation between robustness to input and parameter perturbations within 3-node, 5-node, 15-node, and 30-node TR networks. log 10 G TR I j ð Þ À Á and log 10 G TR p j ð Þ are measures of robustness of a topology j to input and parameter perturbations, respectively, computed from the average over TR networks. Topologies with a low fraction of TR networks (less than 2.3%) are not included. The linear regression for all sizes (3-node, 5-node, 15-node, and 30-node) shows a significant (p,0.0001) correlation between log 10 G TR shows a much steeper slope (1.12 for type2b, 0.33 for type2a, t test = 3.8 and p = 0.0002). This steeper slope may be advantageous for specific biological functions, though both types show strong correlation between the two types of robustness. In the presence of IFF, the two types show no correlation (p = 0.14 and 0.20 for type2a and type2b respectively).

Fine-Grained Analysis within 3-Node Topologies
In this section, we answer the following questions: (1) What is the reason for the large variation around the regression lines in Figs To answer the first question, we speculated that since clearly each of the parameters in a topology will have different robustness values, we might be able to separate the parameters into different categories such that the regression along each category leads to different slope values. If we show this to be true, then as the number of possible categories increases, one expects larger variation in the value of log 10 G TR p for a given log 10 G TR I À Á value. If, in addition, the number of categories is proportional to the number of nodes, then the observed variation would increase for a bigger network, as evident in Fig. 3. In what follows, we investigate this possibility within 3-node networks. Consistent with the 5-, 10-15-, 30-node analysis above, we remove topologies with a low fraction of TR networks (F succ v2:3%) and are left with 2534 topologies to work with.. Next, we separate the parameters of each topology into 7 categories (Fig. 7). Parameters belonging to categories 1 or 2 are those associated with links affecting (i.e., directed towards) the input receiving node, node 1. Those belonging to categories 3 or 5 are associated with links affecting the buffer node, node 2. The rest (in categories 4, 6, 7) are associated with links affecting the output node, node 3. Then, for each category j of a network j,P I , we evaluate the value E p j,P I ,j which takes into consideration only robustness to perturbations in parameters belonging to category j (Eq. 26 in Methods). The corresponding value for the topology j, is G TR p j,j ð Þ (Eq. 28 in Methods). We find that indeed the regression on each of the 7 categories results in a different slope and different correlation strengths. The results of the overall linear regression are shown in Fig. 8A. For the separate categories, we find that the strongest correlation is between log 10 G TR I j ð Þ À Á and robustness to perturbations in parameters belonging to category 1, The second question is related to whether within each topology the parameter subspace corresponding to input robustness is positively correlated with that corresponding to parameter robustness. If they are not correlated, then the two subspaces could be disjoint and the collective/coarse-grained correlation (i.e., the correlation between the log 10 G TR I À Á and log 10 G TR p ) does not support our hypothesis. We follow the same procedure as above and separate the parameters into the 7 categories depicted in Fig. 7. The aim is to be able to compare the results with those in Fig. 8. For each topology j, we perform a linear regression on the relationship between E I j,P  slopes Slope j,j ð Þ, for j~1,2,:::7. The relationship between log 10 G TR p j,j ð Þ and r 2 j,j ð Þ is shown in Fig. 9 while that between log 10 G TR p j,j ð Þ and Slope j,j ð Þ is shown in Fig. 10. As above, the strongest correlations and steepest slopes are found between E I and E p of parameters belonging to category 1, E p j,P I ,1 . For all the topologies, r 2 j,j ð Þ remains $0.9 (Fig. 9A) and Slope j,j ð Þ $0.8 (Fig. 10A). A weaker fine-grained correlation indicates a less collective robustness as indicated by the increase in log 10 G TR p j,j ð Þ (i.e., decrease in parameter robustness) as r 2 j,j ð Þ decreases, for j~1,2,3,5 (Fig. 9A,B,C,E). This pattern does not appear for j~4,6,7 (Fig. 9D,F,G), which is consistent with the results in Fig. 8 Fig. 8H where no correlation is found (as indicated by the high p value). Furthermore, one can map the different clusters appearing in Fig. 8B-H into the clusters that appear in Fig. 10A

Plausible Models of the Role of CheV-P in Bacterial Chemotaxis
The main proteins/receptors involved in E. coli chemotaxis are CheA, CheW, CheB, CheR, CheZ, and CheY. E. coli uses an anticlockwise rotation of its flagella to move forward. A decrease or increase in the concentration of nutrients (chemo-attractants) or harmful chemicals (chemo-repellents), respectively, provokes a change to a clockwise rotation which causes the E. coli to tumble and thus change direction. This signal to the flagella is controlled by the chemotaxis protein CheY. A stimulus (i.e., a change in the chemical concentration in the environment) is sensed by periplasmic binding proteins which couple to CheA in the inner membrane with the help of CheW. An increase in chemoattractant concentrations inhibits the phosphorylation of the receptor complex CheA-CheW (RC-P) (Fig. 11A) while a chemo-repellent enhances it (Fig. 11B). RC-P gives its phosphate group to both CheY and CheB (CheY-P, CheB-P). CheB-P demethylates glutamate residues while CheR enhances methylation. In turn, methylated glutamate (M) enhances the phosphorylation of the receptor complex. The chemotaxis protein CheZ helps speeding the autodephosphorylation of CheY-P [19][20][21] (Fig. 11A-B). For simplicity, we further coarse-grain this network such that M and RC-P interact via a negative feedback loop ( Fig. 11C-D). In the supplementary material (Fig. S3), we demonstrate that there is no significant difference in the results between the topologies shown in Fig. 11A and its coarse-grained equivalent shown in Fig. 11C (slope = 0.79 and 0.77, respectively), though coarse-graining improves the Pearson correlation as it removes the redundant link leading to additional variability. The topology under the influence of a chemo-repellent has a much lower fraction of TR networks and shows no correlation (r = 20.01, p = 0.85) in its un-coarse-grained form (Fig. 11B). It was important to remove the redundancy to obtain a significant correlation (Fig. 11D, slope = 0.37, r = 0.45, p = 10 214 ).
Chemotaxis in many other bacteria is more complex and involves more proteins. One such protein is CheV which generally contains a phosphorylatable domain [24]. Here we consider all possible coarse-grained interactions between phosphorylated CheV (CheV-P), RC-P, and M. The only assumption we make is that RC-P gives its phosphate group to CheV in addition to CheB and CheY (Fig. 11E-F). With this, we obtain 3 3 possible sets of signed directed edges as listed in table 1, where we are considering all 3 possibilities (i.e., activation, deactivation, or no link) for the 3 suggested links.
For each of the 27 topologies, we compute the G TR I and G TR p values (Figs. 12, 13) and the slopes of the regression between E I and E p values of their corresponding TR networks (Figs. 12B,  13B). We compare the results with that of the E. coli topology both under positive (Fig. 12) and negative (Fig. 13) (Fig. 13). Topologies 4,8,18,[25][26][27] are also eliminated as they either show either a negative or no correlation (Fig. 12C-D, 13C-D) under either a chemo-attractant or a chemo-repellent. Topologies 9, 14, 18, 21, 23 are less likely than the rest (13,15,20,22,24) as they have a weaker correlation between E I and E p than E. coli as deduced from the lower p values (Fig. 12C, 13C). Topologies 20 and 22 are less robust to input perturbation than Ecoli when chemo-repellents are the stimulus (Fig. 13), and 24 has a significantly smaller slope. Finally 13 is more robust to parameter perturbations than 15. The distributions of E I , E p for each topology are shown in Figs. S4, S5, S6.

Discussion
In this work, we demonstrated that there is a general positive power-law correlation between environmental and genetic robustness in TR networks, and a statistically significant trend to a directly proportional linear relationship between the two in the limit of large networks. Conversely, monotonically responsive and non-responsive (NP) networks show a weaker relationship than TR ones. Furthermore, this distinction between the two classes becomes more prominent as the size of the networks increases. Therefore, this relationship associated with TR may be relevant to the evolution of biochemical networks. While other factors have played a role in the evolution of genetic robustness, our results show that, for TR networks, as the system evolves to withstand external environmental perturbations, it will, with high probability, concomitantly become robust to certain genetic perturbations.
We speculated that the inverse of the slope is proportional to N {3=4 where N is the number of nodes. We performed the corresponding regression and obtained 1=slope N ð Þ~aN {3=4 zb for a~1:53+0:31 and b~1:01+0:08 (r = 0.9439 , p = 0.008 (1tailed), p = 0.016 (2-tailed)). To confirm our results, we performed a Bayesian analysis for the model slope N ð Þ~1= 1zkN {a with a uninformative flat prior on the parameters and obtained k~1:73+0:31 and a~0:731+0:087 from the second moments of the posterior. Thus the Bayesian analysis confirms the linear regression. For NP networks, the same regression gives  2-tailed)). The value of slope N ð Þ for TR networks in the limit of N large is thus 0.9960.08 while that of NP networks is 0.6860.21. While the latter's regression is not significant at the p = 0.05 level, the value of the intercept did not significantly change for different power values (we tried N {0:64 and N {0:48 ). The statistically significant regression for TR networks implies that as a network evolves to be more robust to input perturbations it will also evolve to be robust to parameter perturbation (and vice versa) at a faster rate. Most importantly, as the size of TR networks becomes larger, the linear relationship between the logarithms quantifying robustness to input and that to parameter perturbations implies that for larger TR networks, G TR P j ð Þis, with statistical significance, and within the computed uncertainty, proportional to G TR I j ð Þ while for larger NP networks, G NP P j ð Þ tends to be proportional to G NP I j ð Þ À Á 0:68 . As standard in the analysis of power-law relationships, we computed the regression using logarithms. For specific biological situations, it may be conceptually more appropriate to compute a direct fit, but for general random networks, we know of no such principle. An exponential fit between G TR P j ð Þ and G TR I j ð Þ for different numbers Figure 10. The slopes within the networks of each 3-node topology divided into 7 categories. Within each topology j, the overall robustness to parameters in category j~1,2, Á Á Á 7 is shown versus slope j,j ð Þ, the slope of the regression line between log 10 E I j,P P and log 10 E p j,P P,j , for j = 1 (A), 2 (B), 3 (C), 4 (D), 5 (E), 6 (F), and 7 (G). of nodes would be difficult to interpret as the power-law is changing with the number of nodes, tending to a constant only as the number of nodes becomes large. A drawback of our method is that the random generation of large networks does not account for reducible topologies which can introduce more variability and thus more error and a lower correlation between the two robustness measures. This makes a comparison between the correlation coefficients of topologies of different sizes a trifle problematic. However, the space of topologies grows so rapidly with the number of nodes that the likelihood of randomly selecting a reducible network decreases precipitously. Similarly, the averaging method does not distinguish between links contributing to the robustness of either input or parameters and those that do not. A method that could pinpoint such links would be useful in this context.
Our results on the adaptability of 3-node motifs differ somewhat from Ref [30] due to our use of a qualitative test, the Pearson shape correlation, for assessing the transient response property of a network. We are not aware of a biologically plausible rationale for an explicit cutoff on the size or speed of a response as biological examples can exhibit both extremes of size or duration of transients. The general motifs shown in the literature [30] need further qualification to be deemed biochemically adaptable. For example, many topologies containing NFL are nonresponsive. Conversely, we show that a subcategory of NFL, NFL1 type2b is particularly robust and exhibits a strong correlation between robustness to input and parameter perturbations (Fig. 6D).
Our results are consistent with biological networks described in the literature. For example, we show that the coarse-grained network topology of E. coli chemotaxis, as described in the literature [17][18][19][20][21], is NFL1 type2b (Fig. S7C), as follows. When the receptor complex is activated, it causes the phosphorylation of the response regulator CheY leading to increased probability of tumbling. An increase in the chemo-attractant level (I) suppresses the activity of the complex and, in turn, the phosphorylation of CheY (Fig. S7A). If I is the input (which we set to always activate the input-receiving node in our computations, for consistency), then we can define the concentration of the input-receiving node as that of the deactivated complex, X 1 (i.e., the activated complex represent X 1 in its deactivated form). In this case, X 1 deactivates CheB which inhibits methylation (M). M activates the complex which is equivalent to deactivating X 1 . The latter inhibits the phosphorylation of CheY (the output) and thus decreases the probability of tumbling (Fig. S7B).
An example of IFF is the Ras model of MAPK cascades discussed in Ref [53]. The input simultaneously activates two factors, SOS and RasGAP which activate and deactivate Ras, respectively and simultaneously. The model is shown [53] to be responsive only when the activation of SOS is faster than that of RasGAP. Thus, one can further coarse-grain it by removing the intermediate node between Ras and the input node (Fig. S8). This reduces to an IFF1 topology.
In Ref [30], all NFL topologies wherein the output node directly affects the input receiving node were found to be not robust or transiently responsive. While consistent with our results showing that NFL7-NFL10 are not robust, note that when the negative feedback loop has a direct and an indirect path, the outgoing and incoming links of the input receiving node must have the same sign for adaptability and parameter robustness to be achieved (see NFL4 and NFL6 as opposed to NFL3 and NFL5 in Fig. 6B). Our work goes beyond pointing out general motifs. We refine subcategories within these motifs and show that, in fact, they do vary in their biochemical adaptation properties.
Traditionally, network motifs represent subgraph topologies that appear in biological networks much more often than one would expect in a randomly constructed network [49], and specific functions were assigned to different types of motifs [41,[46][47][48]. The validity of this approach has been questioned as the frequency of occurrence of these motifs was not statistically significant when compared with corresponding (i.e. same degree) randomly constructed networks [54]. It was argued that one cannot analyze subgraphs independently of the rest of the network as interactions will drastically change the functions assigned to the particular topology [55]. In our work, a motif does not represent a subgraph, rather the topology of the backbone of (possibly much) bigger networks. We use our approach to differentiate between plausible models of the role of the CheV-P protein in bacterial chemotaxis. We find that there are only a few possible ways that CheV-P can be linked to RC-P and M. We suggest that while there are at most 9 possible topologies, the most plausible one has M enhancing the phosphorylation of both CheV and the receptor complex.
Some specific network features have been associated with robustness to environmental variation in bacterial gene expression. Insulating gene expression by different modes of control, from activation to repression depending on the required high or low activity, has been suggested as a general control feature [56].
Our approach to motif discovery can be extended to networks with backbones with more than 3 nodes. While exhaustive enumeration of small motifs with desired functions is fascinating [30,[41][42][43], it is neither immediately evident nor has it been demonstrated in any context that such motifs could be put together to make systems with multiple functions while preserving the robustness or responsiveness properties of the separate motifs. To get to the point where we can plausibly discuss architectural principles in biology, it seems necessary to find general characteristics of classes of networks of all sizes that could perform functions of biological interest. Our work is a step towards this goal.

Notations
Following the same initial setup as in Ref [30], a biochemical network is represented with a directed signed graph S wherein the nodes of the network represent the enzymes. The latter can either be active or inactive and are able to interconvert between the two states. Thus, the elements of the corresponding adjacency matrix A S ð Þ can take the values a ij~{ 1,0,or1, implying that node i deactivates node j, has no effect on j, or activates j, respectively. No parallel links going in the same direction are allowed, i.e., a ij cannot be .1. We divide the nodes into two types, varying nodes and fixed nodes. The latter correspond to inputs and basal enzymes which are added to each network to ensure that each node has at least one activating and one deactivating link. Thus, for an N-node network with N I inputs and N E basal enzymes, A S ð Þ is an n|n matrix where n~NzN I zN E . These concentration values are represented by an n{dimensional vector Assuming that the enzymes are non-cooperative and hence that they obey the Michaelis-Menten kinetics, the rate equations governing the dynamics of the network take the following compact form where h is a unit step function defined as k ij and k m ij are the catalytic and Michaelis-Menten rate constants for the regulation of enzyme j by enzyme i, for j~1, . . . N and i~1, . . . n.
In equation (3), the total concentration of each enzyme is kept constant and normalized (i.e., the concentration of the active form of an enzyme plus that of its inactive form is always equal to one). Thus, 0ƒX i ƒ1 for i~1, . . . N. For all simulations presented here we use only one input, X I 1~0 :5. This particular choice of input concentration should not have a significant effect on our qualitative results, as we have checked explicitly while formulating our hypothesis. Networks are allowed to reach steady state before the concentration of the input is perturbed. We are only concerned with the relative change in steady state concentrations.

Experimental Setup
N-node networks are identified with directed signed graphs S j representing their topology and a set of parameters whereÑ N t is the total number of sampled topologies, excluding those wherein one or more nodes have a total degree of zero or the output node cannot be reached from the input receiving node (Fig. 1). Each topology, in turn, is sampled over a large number of random networks, i.e., a large number of randomly chosen parameter sets P increases exponentially with the size of the networks. For example, N p values for N~3 are around 20, around 40 for N~5, and 500 for N~30. Similarly, the number of iterations (i.e., number of sampled networksÑ N s j ð Þ) required (see Selection Criterion below) also increases with N. For example, forN~3Ñ N s j ð Þ values range between 10 4 and 10 5 , while for N~30,Ñ N s j ð Þ values range between 10 6 and 10 7 . Typically, an iteration takes less than 10 23 seconds of CPU time for small networks (N~3), thus 1 to 2 minutes to test each topology. On the other hand, for large networks (N~30), an iteration typically takes around 0.04 seconds of CPU time, 3 to 5 days for testing each topology.

Pearson Test
We define a TR network as one whose output dynamics has a non-monotonic transient between two steady states as a response to input change (i.e., the steady state values before input perturbation and that after input perturbation). We find the transition time (i.e., the time at which the concentration is maximal/minimal before it starts decreasing/increasing again) and enzyme concentrations, t m and X where X Ã k is the concentration of node k at steady state. We use a Pearson test to determine if a given network is TR. First, we define two functions F ad t ð Þ and F nad t ð Þ as model functions of perfect adaptability and non-adaptability (a monotonically changing network), respectively (Fig. 2): Define corresponding Pearson shape correlations r ad and r nad as where X X N , F F ad and F F nad are the mean values of X N , F ad and F nad . With this, a network is deemed TR if r ad j jw r nad j j. Comparing the absolute values or r ad and r nad instead of the actual values is necessary. Even though our definitions of F ad and F nad will most likely lead to positive values, this is not always the case. The reason is that eq. (6) and (7) assume the perfect case where the differences in the concentrations from the initial steady state have always the same sign (as in Fig. 2 E and F). If instead the difference in concentrations at the transition point is smaller than that at the final steady state (i.e., post-perturbation steady state), then r ad and/or r nad will have negative values. However, that does not matter since we are only interested in the shape of the time-course (see Fig. 2G, for example).
Note that the mean values are taken as the average over all the discretized time-steps; for example, X X N~1 n t X nt t~0 X N t ð Þ for n t time steps. The size of the time-step, dt, is the same for all networks (dt~0:1), but this is not the case for the number of time steps, n t , as the length of the time-course of each network varies depending on how long the network needs to reach a new steady-state (i.e., the rate equations in eq. (3) for all nodes reach zero again after input perturbation. Of course, computationally, the run will stop when the rate equation for all nodes is less than 10 210 ). For example, in Fig. 14A and 14B we show the time-courses (in blue) of two different networks. The network in Fig. 14A needed around 130 seconds (n t~1 300) to reach a steady state, while that shown in Fig. 14B needed around 150 seconds (n t~1 500). Simulations that take too long to reach a steady state (n t w2000) are thrown away and not considered in the analysis (i.e., are thrown away without performing the Pearson test). This cutoff on the maximal number of time-steps allowed is chosen for computational efficiency. Preliminary results showed that for most networks, if a steady state was not reached within 2000 time-steps, it is unlikely it will be reached for a long time. Since we are only interested in the statistical results and since networks are chosen randomly, there is no reason to insist on including a network that takes a lot computational time to reach a steady state. We chose dt~0:1 because we were looking for the largest time-step (to improve computational time) that does not change the statistical results. In preliminary runs, we compared the results of 3-node networks when using dt~0:01 and dt~0:1. The finer time-step allowed more topologies to pass as TR. However, these topologies had very low fraction of TR networks and were removed after the cutoff. Moreover, the statistical results were the same both before and after the cutoff. As mentioned in Experimental Setup above, large networks (30-node) take 3 to 5 days of CPU time for each topology. Using dt~0:01 would increase this simulation time to over a month for each topology which is impractical.

Robustness of the Pearson Test
We test the robustness of the Pearson test described above by comparing the results from the 3-node simulations to those employing instead the Spearman correlation using the same definition of F ad t ð Þ and F nad t ð Þ (Fig. S9A). Both are also compared to simulations using a different definition,F 0 ad t ð Þ and F 0 nad t ð Þ as follows: where a is chosen here to be a~0:3. This new definition allows Environmental & Genetic Robustness Are Correlated F 0 ad and F 0 nad to get to the transition concentration, X m N , at a slower rate, then after the transition point, t m , F 0 nad coincides with F nad while F 0 ad relaxes back to the pre-perturbation steady state, X Ã N at a much slower rate (Fig. 14A,B). In all cases we find no significant difference between the results for 3-node simulations (Fig. S9). This does not mean that there are no variations within individual networks. For example, in Fig. 14, we show the r ad and r nad values of the Pearson test for all the networks corresponding to a typical 3-node topology using both definitions (Fig. 14C). We find that there are 90 out of 21579 networks that were deemed TR in one but NP in the other. A typical time-course where the outcome of the Pearson tests differ or agree are shown in Fig. 14A and Fig. 14B, respectively. In general, most networks do not fall into this category where the values of r ad and r nad are very close such that different definitions of F ad and F nad lead to different outcomes. In Fig. S9, we verify that this change does not affect any statistical observations.

Quantitative Measure of Network Robustness
To quantify the degree of robustness to input and parameter perturbations of a particular network, we calculate the relative change in the steady state concentrations of the output node due to perturbing the input and parameter values, respectively. Let E I and E p be the average of the sensitivity of the steady state concentration of the output to each input and each parameter, respectively. Then, E p j,X I I ,P where the N th node is the output node, P

Quantitative Measures of Robustness of TR Topologies
A robust topology is one that gives rise to robust networks with a higher probability when tested with a large number of parameter sets. Quantitatively, the degree of robustness to input perturbations of a given topology is taken to be the geometric We choose the geometric average as more suitable than the arithmetic average as a conservative approach to detecting a possible correlation, as the latter gives too much weight to much larger outliers.

Selection Criterion
A trial is rejected if it takes too long to reach equilibrium, or its corresponding Jacobian with respect to the node concentrations is singular (i.e., J I I x is noninvertible). With this, we obtain N t |N s matrices for the relative errors E I j,X I I ,P We reject a topology if N s j ð Þ is obviously too small to be statistically significant or F succ j ð Þ takes too long to reach equilibrium (see below), indicating that the parameter space leading to TR networks for that topology is too small.
We sample over 50,000 different topologies for each Nw3 and all possible 3-node topologies (19683), and for each we randomly sample over a large number of parameter sets from a uniform distribution within the ranges k ij [ 0:1{10 ½ and k m ij [ 0:001{100 ½ (whenever a link exists between vertices i and j). For Nw3 the topologies were randomly generated such that the value of each a ij in the corresponding adjacency matrix can take the values 21 or 1 with probability p = 2 each, and a value 0 with probability 1{p. We generated different set of topologies with p~0:1,0:15,0:2,0:3,0:4,0:5. We found no significant difference in the distributions of G I and G p values depending on p. The results shown here represent the collection of all the sets.
We automatically reject trials wherein X m N X I I ?X I I zDX We investigate the effect of the choice of the ranges above by running two separate 3-node simulation. In the first, the parameters are chosen from a uniform distribution in the ranges k ij [ 0:001{2 ½ and k m ij [ 0:001{100 ½ , while in the second, the parameters are chosen from a uniform distribution in the ranges k ij [ 5:0{15 ½ and k m ij [ 0:001{100 ½ . For both ranges we find a significant correlation between robustness to input and parameter perturbations (Fig. S10). For range1 and range2 we obtain the respective values 0.38 and 0.54 for the slopes and 0.73 and 0.68 for the Pearson correlation (Fig. S10A). Furthermore the difference in the slopes becomes insignificant when only networks appearing in both ranges are taken into consideration (Fig. S10B).

Sensitivity Analysis
As discussed above, the degree of input and parameter robustness is seen as inversely proportional to the average of the sensitivity of the steady state concentration of the output to each input and each parameter, respectively. Then, Inserting (22) in (23), we obtain E I j,X I I ,P Similarly, for parameter perturbations, LX Ã N LP a~V P Na and E p j,X I I ,P Fine-Grained Analysis  Figure S2 Correlation between robustness to input and parameter perturbations within TR and NP networks before the cutoff. Correlation between log 10 G I j ð Þ ð Þ and log 10 G p j ð Þ À Á for all topologies showing any number of TR network. log 10 G I j ð Þ ð Þ and log 10 G p j ð Þ À Á are either computed from the average over TR networks (blue) or from the average over NP networks (green). The linear regression for all sizes (3node, 5-node, 10-node, 15-node, and 30-node) shows a significant (p,0.0001) correlation between log 10 G I ð Þ and log 10 G p À Á .   Figure S3 Distribution of E I , E p in the original and coarse-grained Ecoli topologies. (A) Ecoli 1 is the topology shown in Fig.11A and Ecoli 2 is its coarse-grained equivalent shown in Fig. 11C. Their corresponding slopes are 0.79 (r = 0.73) and 0.77 (r = 0.86) respectively. (B) Ecoli 1 is the topology shown in Fig. 11B and Ecoli 2 is its coarse-grained equivalent shown in Fig. 11D. Their corresponding slopes are 20.98 (r = 20.01, p = 0.85) and 0.37 (r = 0.45, p = 10 214 ) respectively. Here we see more variation in the slope than in (A) as the fraction of TR networks is too low for accurate results.
(TIF) Figure S4 Distribution of E I , E p for topologies number 1-3, 5-7, 10-12, 16, 19 when the input is a chemoattractant. These are the topologies that show no TR networks within the sampled parameter space when the input is a chemorepellent. The corresponding slopes, r, and P_values are shown in  Figure S14 The slopes within the networks of each 3node topology divided into 7 categories. Within each topology j, the overall robustness to input perturbations is shown versus slope j,j ð Þ for j = 1 (A), 2 (B), 3 (C), 4 (D), 5 (E), 6 (F), and 7 (G).