Tension and Robustness in Multitasking Cellular Networks

Cellular networks multitask by exhibiting distinct, context-dependent dynamics. However, network states (parameters) that generate a particular dynamic are often sub-optimal for others, defining a source of “tension” between them. Though multitasking is pervasive, it is not clear where tension arises, what consequences it has, and how it is resolved. We developed a generic computational framework to examine the source and consequences of tension between pairs of dynamics exhibited by the well-studied RB-E2F switch regulating cell cycle entry. We found that tension arose from task-dependent shifts in parameters associated with network modules. Although parameter sets common to distinct dynamics did exist, tension reduced both their accessibility and resilience to perturbation, indicating a trade-off between “one-size-fits-all” solutions and robustness. With high tension, robustness can be preserved by dynamic shifting of modules, enabling the network to toggle between tasks, and by increasing network complexity, in this case by gene duplication. We propose that tension is a general constraint on the architecture and operation of multitasking biological networks. To this end, our work provides a framework to quantify the extent of tension between any network dynamics and how it affects network robustness. Such analysis would suggest new ways to interfere with network elements to elucidate the design principles of cellular networks.


Introduction
Decades of experimental studies have established detailed ''wiring diagrams'' of diverse cellular networks. A striking property of many networks is multitasking -the ability to generate different dynamics according to their operating context ( Figure 1A). For example, the mitogen-activated protein kinase (MAPK) pathway involving RAF-MEK-ERK responds to epidermal growth factor (EGF) by triggering transient ERK activation in a graded fashion whereas nerve growth factor (NGF) induced sustained ERK in bistable manner [1]. These tasks directly underlie contrasting biological outcomes: EGF induces proliferation whereas NFG induces differentiation into neurons. Another example concerns the p53 stress response network that mediates arrest, death, and DNA repair functions [2]. In response to ionizing radiation, the network generates multiple pulses of p53 with constant amplitude (i.e. digital) [3] whereas UV-radiation generates a single, broad pulse whose amplitude follows a graded dose-response (i.e. analog) [4]. Insight into how distinct p53 tasks translate into biological outcomes is just beginning to emerge [5].
Multitasking networks are speculated to have arisen through successive elaboration on pre-existing ''core'' processes, representing an evolutionarily feasible route to generate novel biological attributes [6]. Intuitively, reusing a common set of components to multitask can be an economical way to accomplish multiple biological goals. Yet, such a strategy can pose an operational challenge: A dynamic may require network states (each being defined by a set of parameter values) that are ill suited for other dynamics. This concept is related to applications of multi-objective optimization (MOO) algorithms in engineering [7], where two or more, possibly conflicting design aspects are considered. Recently, these approaches have been adopted for biology in problems involving classification, system optimization, and gene regulatory network inference [8]. Here, we use ''tension'' to describe the difference in parameter spaces for distinct dynamics. Intuitively, tension increases with the number of tasks that a network is charged with as each task invariably requires a different subset of parameter values. In the extreme, tension can constrain a network to the point that few additional changes to the network can be tolerated.
A full understanding of network design principles requires an appreciation of where such tensions can arise within networks, their consequences on the robustness of each dynamic, and the strategies used to overcome them. Thus far, however, such concepts have been neglected in quantitative analysis of natural and synthetic pathways. To this end, we have developed a generic computational framework to allow streamlined examination of these questions. We illustrate the use of this framework by analyzing a well-studied RB-E2F network, which plays a pivotal role in regulating cell cycle entry.

The multitasking RB-E2F network
The RB-E2F network has been examined in detail under both normal [9] and pathological [10] circumstances ( Figure 1B and supporting text (Text S1)). In quiescent cells, E2F is silenced by RB [9] whereas E2F expression and activity is modulated by growth stimulation through four ''modules'' -interconnected subsets of the network with a distinct regulatory effect. A sensor module links extracellular growth stimulation and E2F activity: Upon growth stimulation, the MYC level increases and facilitates E2F expression [11] directly and via D-type Cyclins (CYC D ), which potentiate kinases (CDK 4/6 ) to inactivate RB. A positive feedback module (PFB) reinforces E2F expression by two routes: E2F can bind to its own promoter and maintain an activated state [12] and E2F increases expression of Cyclin E (CYC E ) [13], which activates another RB-kinase (CDK 2 ). A negative feedback module (NFB) consists of E2F-regulated genes that include Cyclin A [14] and SKP2 [15] that inactivate E2F binding and induce proteolysis, respectively. Finally, a repression module (R) consists of MYC-regulated genes that down-regulate E2F expression, which may include microRNAs within the miR-17-92 cluster [16] and the ARF tumor suppressor [17].
Three distinct E2F dynamics underlie the response to growth stimuli, depending on the operating context of the network. First, E2F is bistable with respect to serum. Once activated, E2F remains ON even if serum is reduced below the threshold required to activate E2F [18]. In particular, the serum response of E2F exhibits hysteresis, whereby activation of E2F from the OFF state (by increasing serum) and shutting-OFF from the ON state (by decreasing serum) follow different trajectories ( Figure 1B). This property provides a mechanism for cells to enforce two distinct states, quiescence and proliferation [19]: Cells will commit to the cell cycle when a growth stimulus exceeds an activation threshold and to quiescence when signals drop below a maintenance threshold.
Second, E2F exhibits biphasic response to direct MYC stimulation: E2F expression increases with the MYC level when the latter is low, but is repressed when the MYC level is too high [20]. This response restricts the range of MYC levels that can activate E2F. It may represent a safeguard mechanism that allows cells to distinguish physiological levels of MYC induced by serum from transient, potentially oncogenic levels resulting from gene mutation or stochastic gene expression.
Third, in normal cells strongly stimulated by serum, E2F expression exhibits temporal adaptation: It increases to a high level leading up to the end of G 1 before being down-regulated as cells enter the S-phase [21]. As E2F controls expression of many genes involved in DNA synthesis [22], adaptive E2F can both promote coherent induction of DNA replication activities and restrict them to a brief period in S-phase. Indeed, precocious or prolonged E2F activity has been shown to cause replicative stress resulting from deregulated DNA synthesis followed by a DNA damage checkpoint [23,24].

Modeling framework
The starkly different dynamics generated by the same network led us to hypothesize the existence of conflicts that constrain its operation. To examine this issue, we probed several questions by modeling: How (dis)similar are the solution set of parameters that underlie different dynamics? What is the relative difficulty in identifying such parameter sets and what properties do they demonstrate in terms of network performance? In short, for a specific set of dynamics, what is the relationship between tension and robustness?
Here, we have developed a generic computational approach to examine these questions ( Figure 1C). Candidate parameter sets were used to simulate from the model and assigned a score based upon an objective function ( Figure S1A). In a single iteration of the algorithm, randomly initialized parameter sets were subjected to successive rounds of 'mutation' followed by scoring. If a solution was identified, the iteration was terminated or it was terminated without a solution after a defined number of consecutive mutations (in this case 100) without an improvement in the objective score.
This analysis allowed us to enumerate parameter sets that satisfy each particular task (i.e. single) or biologically relevant pairings (i.e. dual). For two tasks (e.g., A and B), tension is calculated as the weighted sum of the log-ratio of median parameter values ( Figure 1D). In the case that each parameter receives equal weighting (i.e. 1/n, where n is the number of free parameters), tension is the average extent each parameter shifts between single tasks. We evaluate robustness according to the ''accessibility'' of dual solutions and ''resilience'' of single-task or dual solutions to parameter perturbation. Accessibility is defined as the fraction of single-task solutions identified as dual. A decrease in accessibility indicates increasing difficulty in locating dual solutions. Resilience is defined as the ability of a solution to maintain some minimal performance after a perturbation (in this case at least 10% of the objective score). This framework can be applied to any kinetic model of cellular networks where objective functions can be quantified.

Tension and coordination between bistable and biphasic responses
We first compared the bistable response to serum and the biphasic response to MYC. From 10,000 iterations we identified a large fraction of solutions for each single task ( Figure 2A). However, only 146 dual solutions were present amongst 4,541 for hysteresis and 14 dual solutions were present in the 4,878 for biphasic. This result corresponds to a dual-solution accessibility of A HB = 0.017, that is, dual solutions represent 1.7% of the total. The rate of solutions identified per iteration and dual accessibility was similar even when only 500 iterations were performed ( Figure S1B), indicating that the result from 10,000 iterations is representative.
Reduced accessibility may reflect tension in the network that arises because single dynamics may adopt disparate states. To examine the correlation between shifts in dynamics and corresponding changes in parameters, we determined the median value of each parameter from all the solutions. By using the values for

Author Summary
Multitasking pervades our daily lives. For example, the technological devices that we increasingly rely upon are now engineered with such multifunctionality or ''integration'' in mind. Similarly, cellular networks also multitask in that they generate multiple, distinct dynamics according to their operating context. Here we show that differences in parameter spaces that underlie different dynamics thus cause a ''tension'', which ultimately constrains network operation. In particular, our analysis reveals that tension negatively impacts robustness by reducing accessibility of parameters able to accomplish two tasks and reduces their ability to withstand perturbations. The presence of tension and its negative impact on network robustness represents a fundamental, generic constraint on the operation of different multitasking networks. The RB-E2F switch, partitioned into four modules: Sensor (cyan edges); repression (R, red); negative feedback (NFB, purple); and positive feedback (PFB, green). This pathway performs at least three dynamic tasks in hysteresis as a reference, we isolated changes specifically associated with biphasic response. According to their influence on each module (i.e., synthesis rates are proportional whereas degradation constants are inversely proportional to module strength), parameters were grouped into four modules (sensor, NFB, PFB, and R). This analysis identified biases in the solution parameters associated with sensor, NFB, and R modules ( Figure 2B), whereas changes to the PFB parameters were divergent (see supporting Text S1, Discussion). Note that the overall distribution of NFB values were quite similar between different dynamics despite a change in median ( Figure S2A). The changes across all median parameter values resulted in a tension of 0.34 (i.e. average shift in parameter value) between hysteretic and biphasic tasks.
Parameter distributions may be highly irregular, raising the issue of how the median may perform as a summary of each solution set. An alternative approach to compare distributions is to calculate the Kullback-Leibler (KL) divergence (Text S1). Consistent with the results obtained using median values, the largest KL divergence involved parameters of R ( Figure S2B). More subtle distances in NFB, Sensor, and NFB were also present. In this case, the tension value (0.08) was calculated as the average KL divergence. All subsequent analyses were done by using the median values.
The difference between the two dynamics can be largely accounted for by the strength of sensor and R modules -the product of free parameters constituting each module ( Figure 2C). Given this observation, an effective strategy to reconcile the tension is to dynamically configure these modules: Increasing their strengths would favor biphasic response, while decreasing them would favor hysteresis. In contrast, changes in other modules would be less critical. We term this dynamical 'network reconfiguration'.
The overlap between R and sensor ( Figure 2C), however, also suggests the possibility to accommodate the two dynamics by using common parameter sets, which, by definition, represent dual solutions. We performed 10,000 search iterations using a composite objective function that represents the product of hysteretic and biphasic objectives (see Text S1, Materials and Methods) which allowed us to identify an additional 1,290 dual solutions ( Figure 2A). Most dual solutions were concentrated in the overlap between single-solution sets, consistent with the notion that they represent a hybrid of parameters from single dynamics ( Figure S2C). To validate the distribution of these dual solutions, we also attempted a search using a dual objective function composed of the sum of individual objectives. In addition, we performed a search with single hysteretic and biphasic solutions as a starting point, mimicking the successive elaboration of network tasks. In each case, the distribution of solutions parameters was virtually indistinguishable ( Figure S2 A and C). This supports the notion that the distribution of dual solutions is representative.
Simulations show that a typical dual solution could indeed generate both dynamics ( Figure 2D). Consistent with Figure 2C, weakening the R module (by substituting it with the median value from hysteretic solutions) diminished the repression of E2F at high MYC, thus diminishing the biphasic response. In contrast, strengthening the R module (by substituting it with the median value from biphasic solutions) maintained the biphasic response to MYC but eliminated hysteresis by weakening overall E2F response. Weakening the sensor shifted the hysteretic response to higher serum inputs but diminished the E2F levels achieved in response to MYC ( Figure S2D); strengthening the sensor eliminated hysteresis and broadened the biphasic response by stimulating an increase in E2F at relatively low doses of input.
A caveat of such dual solutions is their reduced accessibility ( Figure 2A). In addition, it is interesting to examine if tension could also impact their resiliency to perturbation. To examine this, we selected fifty representative solutions from each category in the vicinity of their respective medians, subjected each one to 10,000 parameter perturbations, and determined the fraction of perturbations that retained at least 10% of the initial objective score. This analysis revealed that biphasic response was a more resilient property than hysteresis overall ( Figure 2E and Figure S2E).
Although the median resiliency of dual solutions was slightly lower than single solutions, this change was not significant, suggesting this tension had a minimal impact on the performance of dual solutions. As such, properly configured sensor and R modules can accommodate both dynamics. This could be achieved by engaging the R module only when MYC is sufficiently high, yet simultaneously enhancing the sensitivity of E2F to MYC stimulation. This notion is consistent with the distinct modes of MYC regulation in physiological and pathological contexts. Physiological stimulation, e.g., by serum, of arrested cells leads to a pulse of MYC that drops to a low level throughout the cell cycle [11], which is unable to trigger the R module. Still, a strong sensor module would enable robust generation of E2F switching behavior despite relatively low MYC levels (second column of Figure 2D). In contrast, more elevated and persistent levels of MYC, due to overexpression or stochastic gene expression, would trigger the R module and result in biphasic response.

Tension and coordination between hysteretic and adaptive responses
Using the same approach, we found that the accessibility of dual solutions involving hysteretic and adaptive dynamics was 7-fold lower compared to biphasic behavior (A HA = 0.0024 compared to A HB = 0.0170) ( Figure 3A). This decrease was accompanied by an elevated tension between hysteresis and adaptation (T HA = 0.48 compared to T HB = 0.34). Compared to hysteresis, adaptation is associated with parameters defining moderately enhanced sensor and R modules, and a drastically stronger NFB module ( Figure 3B and Figure S3A). Changes in parameters associated with PFB were without coherent bias (Text S1, Discussion). Consistent with these results, the dominant shift in KL divergence involved NFB parameters ( Figure S3B). Furthermore, the tension (average KL divergence) between hysteretic and biphasic dynamics (0.08) is lower than that between hysteretic and adaptive dynamics (0.11). These observations suggest that an effective strategy to reconcile the drastic tension is to dynamically configure these modules, particularly for the NFB: Increasing its strength favors adaptation, while decreasing it favors hysteresis ( Figure 3C). Reflecting their 'hybrid' nature, dual solutions were concentrated in the overlap between individual dynamics when plotted as a function of sensor and NFB strengths ( Figure S3 A and C).
To examine the specific contribution of NFB and sensor in modulating these dynamics, we varied its strength in a typical dual solution. Simulations confirmed its ability to generate both dynamics ( Figure 3D). Weakening the NFB module (by substituting it with the median value from hysteresis solutions) eliminated response to growth signals. (C) Algorithmic approach to search parameter space. (D) Calculating tension. Given a network, solution sets of parameters able to generate each dynamic (A and B) are identified. Tension between tasks (T AB ) is defined as the weighted sum of the log-ratio of median parameter values. doi:10.1371/journal.pcbi.1002491.g001 the adaptive response and strengthening it (by substituting it with the median value from adaptive solutions) increased the precision of adaptation, consistent with its requirement for this behavior [25]. Weakening the NFB module also enhanced hysteresis to the point that E2F expression became irreversible (i.e. the solid and dotted curves do not meet at low serum). Yet, strengthening it diminished hysteresis by interfering with maintenance of the E2F ON state upon reduction in serum. The sensor strength had a more general impact ( Figure S3D). The sensor strength for the dual solution seemed to be near optimal for hysteresis; either weakening or strengthening it led to almost elimination of hysteresis.
The strong tension between dynamics corresponds to a greatly reduced dual accessibility and suggests that they may be operational over a much restricted parameter space ( Figure 3C). Indeed, here tension penalized the performance of individual dual solutions: Dual solutions were significantly less resilient to perturbations than single solutions in maintaining both hysteretic and adaptive dynamics ( Figure 3E and Figure S3E).
The drastically reduced accessibility and robustness of dual solutions suggests that they would be ineffective in accommodating both dynamics. Instead, dynamic network reconfiguration is likely critical, which is consistent with the operation of the network: the negative feedback on E2F has a significant time-delay in its operation. In the G 1 phase, the Anaphase-Promoting Complex/ Cdh1 (APC Cdh1 ) keeps negative feedback from both CYC A [26] and SKP2 [27] low by targeting them for proteasomal-mediated degradation. Upon progression to G 1 /S, E2F activity increases and induces CYC E -which is resistant to APC Cdh1 -engaging sole positive feedback. SKP2 and CYC A levels are eventually allowed to increase through E2F-mediated induction of Emi1 [28] which targets APC Cdh1 for destruction. This is reinforced through positive feedback as CYC A itself can also target APC Cdh1 for destruction [29]. Delay is also achieved at the transcriptional level through ordered release of Cyclin E and Cyclin A from RB-mediated repression [30]. This temporal coordination has been speculated to enforce a brief time window between DNA replication origin licensing mediated by CYC E and origin deactivation and initiation of DNA synthesis mediated by CYC A [31,32]. Sequential triggering of positive and negative feedback appears to be a generic, systems-level organizational principle of networks underlying cell cycle control conserved throughout evolution [33,34]. Our analysis suggests an additional role for the temporal coordination: it represents dynamic network reconfiguration that accommodates robust hysteretic and adaptive E2F responses.
In addition, the tension between dynamics can potentially be alleviated by increasing network complexity. For example, eight E2F members of the E2F family have been identified in mammals; some members can functionally substitute for one another [35]. E2F1 and E2F3 are part of the ''activator'' subgroup required for cell cycle entry of fibroblasts from quiescence [36]. We wondered if such apparent redundancy could reduce tension. To test this notion, we extended our model to include an additional E2F member (E2F9) expressed in parallel with E2F ( Figure 4A and Text S1, Mathematical Model). In particular, the model includes distinct parameters governing production and degradation of each E2F copy. On the other hand, we assumed that the biochemical activity of each E2F copy was indistinguishable and could contribute in an additive manner to overall E2F output (i.e. hysteresis and adaptation) as well as to downstream gene expression (i.e. Cyclin E and Cyclin A) via shared parameters.
The added complexity indeed led to a 3.1-fold increase in dual solution accessibility (A 2xE2F HA = 0.0052 compared to A HA = 0.0024) ( Figure 4B). This was accompanied by a reduction in network tension with dual E2F (T 2xE2F HA = 0.39 versus T HA = 0.48) ( Figure 4C). This is reflected in the more modest extent to which the NFB and R modules shifted between hysteretic and adaptive dynamics ( Figure 4C and Figure S4A) and the greater extent of overlap in their distributions ( Figure S4B). Importantly, inclusion of an additional E2F copy was sufficient to increase the resilience of dual solution adaptation such that the median was not significantly different from single task solutions ( Figure 4D). In contrast, this additional complexity did not have a significant impact on the resilience of hysteresis associated with dual solutions. Why this fragility of hysteretic dynamics persists in such dual solutions is not clear. Nevertheless, these results are consistent with the notion that increasing network complexity reduces tension and the corresponding penalty on some aspects of robustness.

Discussion
Quantitative modeling has been widely adopted to examine design principles of biological networks. Many studies have provided important insight into the ways networks generate particular dynamic responses [25,37]. To date, however, how a multitasking cellular network accommodates different dynamics is poorly understood, despite the recognition of their wide presence and importance. Here we have developed a general approach to quantify tension between different dynamics, which we have applied to a well-established network underlying cell cycle progression. In general, our analysis is consistent with an inverse relationship between tension and the overall robustness of network operation ( Figure 5). In the face of moderate tension, common or 'one-size-fits-all' parameter sets could be attractive as they avoid the need for additional, possibly complex, mechanisms to coordinate system parameters. However, dynamic network reconfiguration may be critical to resolve strong tension. Though dual solutions exist, there is a pronounced penalty on the accessibility and resilience of these solutions. Our approach is general in that it can be applied to any other network with behaviors that are distinct and quantifiable. In the case where a network demonstrates numerous tasks (including the RB-E2F network), accessibility, tension, and resiliency can be reported by an ''adjacency matrix'', reporting all interactions in a pair-wise fashion.
Our findings have several implications for our understanding of the RB-E2F switch as well as a variety of other multitasking networks (Table 1). First, our generic framework provides additional criteria to assess model selection, sometimes favoring choices that are not intuitive. In the case of the RB-E2F network, the relatively high tension between hysteretic and adaptive tasks suggests a critical need for additional mechanisms able to delay replacing it with median value from solutions for hysteresis (2R) or biphasic response (+R). E2F is expressed in mM. (E) Resilience of individual solutions. Boxplots summarize the resilience of 50 solutions for hysteresis, biphasic, and Dual (H/B). Solutions were selected by identifying the smallest box (values of sensor and R) centered on the median containing 50 solutions. Resiliency is defined as the ability to maintain at least 10% of their objective score following parameter perturbation. Y-axis shows the fraction of 10,000 repeated perturbations to a particular solution that are resilient. Circle indicates median; Medians are significantly different at the 5% significance level if there is no overlap between intervals defined by their triangular notches. doi:10.1371/journal.pcbi.1002491.g002 negative feedback (i.e. CYC A ) or buffer parameter changes (i.e. E2F duplication). Another example involves a study by Ashall et al. [38] concerning how different pulsatile TNF-a input patterns encode unique NF-kB nuclear translocation dynamics ( Figure  S5A). The authors were prompted to propose an alternative network model wiring when they were unable to find common parameter sets that could satisfy all NF-kB tasks using a traditional model. Using their data, we calculated that the tension between two tasks (''Continuous'' and ''60 minute'') was reduced from T Con/60 = 1.69 to T9 Con/60 = 0.84 in the alternative model along with a corresponding increase in accessibility from A Con/60 = 0 to A9 Cont/60 = 0.29 ( Figure S5B). What system-wide values of tension and accessibility are across all model parameters remains to be seen. Nonetheless, our study suggests that dynamic shifting of parameters is more desirable from the perspective of robustness.
Second, tension has the potential to affect network evolvability [6]. In particular, coopting additional functions could interfere with pre-existing network dynamics (i.e. partial overlap of solution space), thereby reducing the ability of the network to tolerate additional alterations. For example, Meir et al. [39] modeled the ability of the Notch-Delta signaling network [40] to generate three spatial cell fate patterns -''2-cell'', ''7-cell'' and ''Line'' -attributed to the pathway during animal development. They showed that the solution spaces for these tasks were only partially overlapping ( Figure S5 C and D): Only 25% of solutions for the ''2-cell'' tasks could accommodate a ''7-cell'' pattern while nearly 80% of ''7cell'' solutions could also produce ''2-cell'' patterning corresponding to an accessibility of A 2-7 = 0.51 for dual solutions. Also, parameters for ''Line'' overlap to an even lesser extent with solutions for the two other tasks. From this, the authors speculated that existence of universal parameter sets represent an evolutionarily feasible route towards the goal of achieving novel functions. On the other hand, these same observations offer direct support for our argument that tension reduces robustness and constrains a network's capacity to adopt additional tasks. An intriguing possibility is that dynamic shifting, increased complexity, or other strategies may enable a network such as this to increase its workload [41].
Third, by focusing on the coordination of different tasks, our methodology can provide novel, experimentally testable hypotheses concerning what mechanisms are tied to potential conflicts between dynamics and how they are resolved. For example, Santos et al. [42] showed that the MAPK cascade, consisting of RAF, MEK1/2, and ERK1/2, demonstrates distinct dynamics and contrasting phenotypes in response to EGF and NGF (Figure 1 and Figure S6A). Importantly, EGF stimulated negative feedback between ERK and RAF whereas NGF stimulated positive feedback. The growth-factor context-dependent MAPK topologies are a clear example of tension between dynamics and the functional role of network reconfiguration. Analogous to our results showing that modulating NFB strength could impact hysteresis ( Figure 3D), small-molecules used to constitutively suppress and sustain positive feedback could swap ERK dynamics and physiological effects of NGF and EGF. Furthermore, the authors showed that partial activation of positive feedback via interfering RNA (RNAi) generated an intermediate ability of EGF to induce differentiation, suggesting a quantitative relationship between tension and phenotypic outcome.
Another example involves the multifunctional response of the p53 tumor suppressor. Batchelor et al. [4] demonstrated that repeated, digital pulses stimulated by c-radiation (c-IR) required WIP1-mediated negative feedback whereas UV radiation generated a single, graded p53 response ( Figure S6B). Importantly, suppression of negative feedback by RNAi against WIP1 was sufficient for c-IR to generate a p53 response characteristic of UV [43]. These observations represent a clear demonstration of tension between dynamics attributed to negative feedback, and its reconciliation through duplication and diversification of the network (i.e. ATM and ATR). In retrospect, our framework provides a rational means to identify such network tension, which may not easily arise from intuition alone or even a deep knowledge of the network, especially when tension arises from subtle and/or multiple parameter shifts.
For the RB-E2F network, our detailed examination of tension and robustness provide experimentally testable hypotheses. First, our results suggest that the strength of negative feedback acting upon E2F is inversely related to the extent of hysteresis. The strength and timing of NFB could be realized through a smallmolecule inducible Cyclin A expression construct. Alternatively, premature Cyclin A activity could be achieved through introduction of an N-terminal deletion mutant resistant to APC/Cmediated destruction [44]. The effect of this on the E2F doseresponse to serum could be readily achieved using a previously devised fluorescent reporter for E2f1 [45]. Second, this same experimental system could be used to test the hypothesis that additional copies of E2F insulate the hysteretic response from premature or intensified NFB. Finally, our results lead directly to the hypothesis that strong NFB will reduce the robustness of networks able to accommodate both dynamics. Such a question would be best suited using a synthetic biology approach and predicts that circuits with both bistable and adaptive dynamics would arise with relatively mild NFB.  If tension places a fundamental constraint on the operation and architecture of multifunctional networks, it has implications for engineering of synthetic biological systems. To date, most efforts have focused on engineering of gene circuits with limited, dedicated functions. More complex functions can then be realized by integrating well-defined modules [46,47]. For those functioning in individual cells, however, this strategy is limited by the ability to insulate different modules as well as the inevitable burden they impose upon cells which can undermine desired functionality [48]. As such, it may be more effective to explore strategies that include dynamic network reconfiguration to perform multiple functions in a robust manner. In this case, synthetic biology may take a cue from nature: Rather than attempting to generate an everexpanding toolkit of biological components, an emphasis will be placed back upon the vast potential in differential regulation of existing entities.

Numerical simulations
Simulations were performed with Matlab, version R12 (Mathworks, Natick MA) employing the ode15 solver. Figure S1 Objective functions for search algorithm and convergence. (A) Objective functions used to quantify numerical simulation output. (Left) Hysteresis is defined as a minimal path difference (DP = 0.5) in E2Fm at 24 hours after an increase in serum from 0.01% or decreasing from 10%. This was calculated by applying the Matlab function trapz to the difference in steadystate E2F values generated by decreasing and increasing serum. (Center) The relative adaptation in E2Fp was calculated by DF/DI over 25 hours and a minimum threshold of 0.80 defines a solution. DI is the difference between initial and peak levels and a minimal DI is enforced to filter out trivial solutions. DF is the difference between peak and final levels. (Right) Biphasic behavior is defined by the extent of E2Fm suppression relative to initial increase (DF/ DI) at 36 hours after a change in MYC synthesis rate (parameter ke MYC ). A minimum threshold of DF/DI = 0.80 defines a solution. A minimal absolute value of DI also applies in this case. It should be noted that hysteresis, adaptation, and biphasic behavior could be measured at the protein level without loss of generality.  Figure 2D. Sensor strength was decreased and increased by substituting median value for hysteresis (2Sensor) and biphasic (+Sensor), respectively. The value of d E2Fp 21 was increased (+d E2Fp

21
) by using the median value from hysteresis. (E) Evaluation of resilience for representative solutions. The change in objective score relative to original is plotted as a function of total parameter variation (K). Shown are results of 10,000 perturbations. Resilience of a perturbed parameter set is that maintaining at least 10% of its score (i.e. log value greater than 21).  Figure 3D. Sensor strength was decreased and increased by substituting median value from hysteresis (2Sensor) and adaptation (+Sensor), respectively. The value of k E2Fm and d E2Fp 21 were increased (+k E2Fm ,d E2Fp 21 ) by using the median value from hysteresis. (E) Evaluation of resilience for representative solutions. See legend for Figure S2E  The NF-kB pathway mediates stress signals including those from the TNF-a cytokine. Cells treated with different temporal patterns of TNF-a given in 5 minute pulses display distinct NF-kB dynamic 'tasks'. (B) A previous 'traditional' model of the pathway is unable to accommodate all tasks with a common parameter set whereas an alternative, ''Triple-feedback'' model with Ikk feedback is able to. Tension and accessibility of dual solutions involving the continuous and 60 minute TNF-a pulsing protocols were calculated from data using the A20 degradation rate parameter. Adapted from Ashall et al. [25]. (C and D) Tension in the Notch-Delta multitasking network. (C) Notch-Delta signaling leads to differential expression of Achete (AC)/Schute (SC) and binary cell fate patterning in adjacent cells during fruit fly development. The network is able to translate an initial pattern of AC/SC expressed at moderate levels into a final ON/OFF pattern. (D) Calculation of accessibility of dual 2-and 7-cell pattern. Orange bars indicate subset of parameters for each single task that are dual. Accessibility calculated from data presented by Meir et al. [26]. (TIF) Figure S6 Multifunctional networks with diverse, context-specific dynamics. (A) The MAPK pathway multitasks. Stimulation of neuronal precursor cells with EGF and NGF elicit distinct dynamics and translate into opposite phenotypic out-comes. Protein Kinase C (PKC) is required but not sufficient for positive feedback. Small molecules used to sustain positive feedback (phorbol-12-myristate-13-acetate (PMA)) or preclude it (Go7874) were sufficient to swap EGF-and NGF-mediated dynamics and cellular outcomes. Adapted from [27]. (B) The p53 stress response pathway multitasks. p53 can mediate cell stress signals and controls the expression of genes that mitigate their effects. Double-strand (ds) breaks induced by ionizing radiation induce recurrent pulses of p53 that are whose amplitude is doseindependent [28]; Single-strand (ss) DNA adducts induced by UV cause a large pulse of p53 that is graded in terms of peak response [29]. Colored links indicate interactions (i.e., synthesis rate in blue and NFB in red) activated in a stimulus-specific fashion. The p53 network has also been shown to respond to a large panel of cell stresses and other physiological contexts, with dynamics that are poorly understood. (TIF) Text S1 Detailed description of modeling and mathematical framework. The supporting text opens with a discussion section describing the theoretical model of the RB-E2F switch underlying mammalian cell cycle control along with a discussion of the role of positive feedback module in the adaptive and biphasic E2F responses. Following this is a materials and methods section that describes the computational approach used to identify parameter solutions that satisfy RB-E2F network dynamics. Concluding this material are definitions of tension, Kullback-Leibler divergence, and measures of robustness. (DOC)