Research article | Open | Published:
Constructing network topologies for multiple signal-encoding functions
BMC Systems Biologyvolume 13, Article number: 6 (2019)
Cells use signaling protein networks to sense their environment and mediate specific responses. Information about environmental stress is usually encoded in the dynamics of the signaling molecules, and qualitatively distinct dynamics of the same signaling molecule can lead to dramatically different cell fates. Exploring the design principles of networks with multiple signal-encoding functions is important for understanding how distinct dynamic patterns are shaped and integrated by real cellular networks, and for building cells with targeted sensing–response functions via synthetic biology.
In this paper, we investigate multi-node enzymatic regulatory networks with three signal-encoding functions, i.e., dynamic responses of oscillation, transient activation, and sustained activation upon step stimulation by three different inducers, respectively. Taking into account competition effects of the substrates for the same enzyme in the enzymatic reactions, we searched for robust subnetworks for each signal-encoding function by three-node-network enumeration and then integrated the three subnetworks together via node-merging. The obtained tri-functional networks consisted of four to six nodes, and the core structures of these networks were hybrids of the motifs for the subfunctions.
The simplest but relatively robust tri-functional networks demonstrated that the three functions were compatible within a simple negative feedback loop. Depending on the network structure, the competition effects of the substrates for the same enzyme within the networks could promote or hamper the target functions, and can create implicit functional motifs. Overall, the networks we obtained could in principle be synthesized to construct dynamic control circuits with multiple target functions.
The relationship between structure and function in biological circuits is a focus of systems biology [1, 2]. Structure reflects the topology of interactions within a circuit, and function is usually considered a static quantity (a steady state) or a temporal behavior (dynamics) of the circuit’s output . In various signaling systems, cells transmit information via dynamic control of signaling molecules . The information about the identity and quality of a stimulus can be encoded in temporal patterns of a signaling molecule by modulating its amplitude, duration, frequency and so on. Cells mediate gene expression programs and induce cellular responses according to the sequential dynamics of signaling molecules [3,4,5,6,7,8,9,10,11,12,13]. Furthermore, signaling dynamics usually have greater information transmission capacities compared to non-dynamics responses [14,15,16]. For example, the tumor suppressor p53 shows a series of pulses with number modulation under γ-radiation, resulting in cell-cycle arrest. By contrast, UV radiation triggers a single p53 pulse with a dose-dependent amplitude and duration, leading to apoptosis [3, 17, 18]. In yeast, the transcription factor Msn2 plays a major role in the general stress response. Glucose limitation stress induces nucleocytoplasmic shuttling of Msn2. The duration of the initial peak and frequency of the sequential bursts in Msn2 increase with the strength of the stress [6, 19, 20]. Under osmotic stress, Msn2 exhibits a similarly initial nuclear peak and returns to normal. The duration of the peak increases with the osmotic pressure [3, 21,22,23,24,25]. In addition, oxidative stress leads to prolonged nuclear Msn2 accumulation whose amplitude is dependent on the concentration of H2O2 [6, 26, 27]. Other systems, such as ERK, NF-κB and Crz all drive distinct responses via dynamic control [3, 28,29,30].
These observations in real signaling systems raise fascinating questions about how distinct input information is encoded into different dynamic patterns of the signaling molecule and the design principles for such networks with multiple signal-encoding functions. Instead of considering the complex molecular interaction details in real signaling networks, we here intend to perform a theoretical investigation of multi-node enzymatic regulatory networks with three signal-encoding functions, i.e., dynamic responses of oscillation, transient activation, and sustained activation when a respective input node is stimulated. These dynamic patterns are similar to the observed dynamic response of Msn2 in yeast under glucose limitation, osmotic stress, and oxidative stress, respectively.
General computational search strategies for detecting network topologies for a target function include evolutionary algorithms and the network enumeration approach . These approaches have been used to explore simple network solutions for dynamic functions such as oscillations [31,32,33], switch-like responses [34, 35], adaptation [36,37,38] and dose-response alignment . Core motifs and design principles for these individual dynamic functions have been investigated previously. By contrast, multi-functional networks that can execute several distinct functions have rarely been considered , despite the accumulation of experimental observations of multiple signal-encoding processes. The stumbling blocks for searching multi-functional networks are mainly due to the increased network scale and the functional complexity, which require enormous computing power. As the multi-functional networks we consider have at least four nodes, with three nodes for inputting three external signals and one node for output, it is practically infeasible to enumerate all possibilities to detect multi-functional network topologies. Evolutionary approaches are also inapplicable for deducing the architectural landscape of tri-functional networks. Accordingly, here we adopt the scenario of module combination for searching multi-functional networks capable of the desired multiple signal-encoding function. The module combination approach has recently been used to generate bi-functional networks by hybridizing simpler mono-functional modules . It was also applied successfully to detect network structures with Pavlovian-like functions capable of learning and recalling .
In our calculations, the external stimulus, which is a temporal step function, acts on three different input nodes, and generates on the output node sustained oscillations (stable limit cycle), transient activation (adaptation), and sustained activation, respectively. This mimics the signal-encoding of Msn2 in yeast in response to glucose limitation, osmotic stress, and oxidative stress, respectively. To obtain the target multi-functional networks, we search first for robust subnetworks (sub-functional modules) for each target input-output response by enumerating all three-node network topologies and then meld different modules through node merging. All computations are performed with enzymatic regulatory networks. In addition, we take competition effects among different substrates into account in situations when two or more types of substrates bind competitively to the same enzyme. This ubiquitous hidden interaction has been considered in gene regulatory networks , but is commonly omitted in most theoretical or experimental models in enzymatic circuits . We find that competition effects can increase the diversity of functional modules by creating implicit cross couplings within the circuits, whereas in combined networks, competition effects can promote or hamper the overall functionality. These findings may be helpful for understanding multi-functional networks in real biological systems and also provide a guide to construct information-encoding circuits in synthetic biology.
Module selection for each signal-processing function
We use the method of three-node-network enumeration to explore the functional networks of oscillation (F1) and adaptation (F2), respectively. In a two- or three-node network, one node serves as the input node (Node R) and one as the output node (Node O), and a third node will be the intermediate node (Node M). For each topology, we sample 10,000 parameter sets using the Latin hypercube sampling method . For each parameter set, the ordinary differential equations that describe the network dynamics in response to one of the two stimuli are then solved numerically. The topologies with at least one or more parameter sets that produce the targeted dynamical function are considered functional circuits. The robustness of sub-functions F1 and F2 is measured by the number of parameter sets, i.e., Q1- and Q2-values, that generate the target dynamic functions. For functions F1 and F2, we select the network topologies with relatively high Q-values and fewer edges as functional F1 and F2 module pools.
The oscillation (F1) pool contains 73 networks with 3 to 6 edges with Q1 > 50 and 8 simplest networks with 3 edges with 10 < Q1 < 50 (as marked in Fig. 2a). Clustering analysis of these networks (Fig. 2b) shows that the basic oscillation motifs for achieving function F1 are negative feedback loops between the three nodes. Motifs a1 and a3 are simply two forms of the repressilator, and motifs a2 and a4 are the delayed feedback loop. All 73 networks contain explicitly one of the two negative feedback loops which are well-known motifs for oscillations. This is true even for networks with edge numbers greater than 6. Additionally, the F1 functional pool also contains the 8 simplest networks. They do not explicitly have the structure of a repressilator or delayed feedback loop (Fig. 2c). These topologies achieve self-sustained oscillations because one of the two basic motifs of the repressilator and delayed feedback loop is implicitly contained in their topologies due to the hidden interactions. Two such examples are illustrated in Fig. 2g, where the oscillations are achieved through the implicit a3 and a4 motifs. Thus, the modular pool of F1 consists of 81 networks as marked in Fig. 2a.
Similarly, as marked in Fig. 2d, the F2 functional pool also consists of 81 selected networks, i.e., 76 networks with 3 to 6 edges and Q2 > 55, and 5 networks with 3 edges and 10 < Q2 < 55. Clustering analysis of the 76 networks reveals the backbones of the functional circuits (Fig. 2e), two negative feedback loops with a buffer node (motifs b1, b3), a delayed feedback loop (motif b2), and an incoherent feedforward loop with a proportional node (motif b4). The motifs for function F2 are similar to Ma’s results , but different in that motifs b3 and b4 (for which Q2 < 55) are less robust than motifs b1 and b2 (for which Q2 > 55). This difference is attributable to the inclusion of competitive regulations in both motifs b3 and b4: node M3 and node O3 compete for enzyme R3 to convert from the inactive to active state. Consequently, there exit implicit positive interactions between nodes M3 and O3, and the regulatory interactions are more complex than those discussed in , in which the competitive effect was ommitted. Figure 2f illustrates the 5 networks with three nodes and lower Q2 value as marked in Fig. 2d. These simple functional circuits all have implicit interactions, of which two are motifs b3 and b4. The left three motifs have very different topologies from the ordinary adaptation motifs. Checking the implicit interactions reveals that these seemingly very different adaptation motifs all have implicitly the core structures discussed in detail in . Two such examples are illustrated in Fig. 2h where adaptation is achieved through the core structures of b1 and b3.
As demonstrated in Fig. 1 and Fig. 2, the inclusion of competition effects into the circuit dynamical description generates mutual activation or inhibition interactions. Taking these implicit interactions into account benefits analysis of the functioning of networks when we judge whether a specific dynamical behavior is possibly supported in a network topology. It has been well known that oscillations are generated in circuits with topologies of negative feedback loops and that adaptation is typically achieved with topologies of incoherent feedforward or negative feedback loop with buffering node. For the circuits that do not have apparent topologies that support oscillations or adaptations, one can add the implicit interactions onto the network and then check whether the newly formed network have the core structure of oscillation or adaption. In this way, implicit functional motifs that are previously unknown can be found as depicted in Fig. 2g and Fig. 2h.
The functional circuits with sustained activation (F3) are not searched by enumeration. We prefer to adopt the simplest structure, i.e., a link of activation from the input node R3 to the output node O3. Thus, we have constructed the module pools for the sub-functions of oscillation (F1, with 81 topologies), transient activation (F2, with 81 topologies), and sustained activation (F3, with 1 topology), respectively.
Module combination for bi-functional networks of F1 and F2
By combining sub-function networks in the F1 and F2 module pools, we next construct F1-F2 bi-functional networks that encode the external signal received by nodes R1 and R2 respectively into oscillations and adaptation in the output node O. The network topologies from the two sub-functional pools are united by merging their nodes. As demonstrated in Fig. 3, we adopt the following node-merging rules:
Output O nodes have to be merged with each other as we assume a unique output node.
M nodes may be merged with each other (unless there are contradictory linkages).
R nodes may be merged with M nodes from the other pathway (unless there are contradictory linkages).
R nodes for receiving signals cannot be merged with each other.
R nodes are not permitted to merge with O nodes.
Here we do not consider the merging between R and O nodes and node merging within pathway. As there are two input nodes for receiving the signals and a unique output node in a bi-functional network, there are three possible ways to combine the F1 and F2 circuits as illustrated in the upper left of Fig. 3a-c. One approach is to simply merge the output nodes O1 and O2. Another is two-node merging of M1-M2, R1-M2, and R2-M2 in addition to the O1-O2 merging. The third approach is the three-node merging of O1-O2, R1-M2, and R2-M1.
We enumerate all possible networks combined from the two module pools as candidate topologies for bi-functional circuits. For each candidate network, we randomly sample 10,000 parameter sets and check whether it is bifunctional. The number of parameter sets that can achieve both oscillation and adaptation is recorded as the Q1,2-value to estimate the robustness of the bi-function. In the bi-functional networks with Q1,2 > 0, we obtain 3215 five-node, 1337 four-node, and 41 three-node networks by one-, two-, and three-node merging, respectively. The Q1,2 distributions of these three classes of bi-functional combined networks are illustrated in the lower left of Fig. 3a-c. The networks with relatively high Q1,2 are highlighted with red dashed boxes in the histograms of Fig. 3. These robust networks obtained by one-, two- and three-node merging are then clustered (middle panels in Fig. 3a-c). Structure skeletons from the clustering analyses are listed alongside (right panels in Fig. 3a-c). The robust bi-functional networks obtained via one-node and two-node merging are composed of hybrid motifs from the F1 and F2 module pools, such as hybrids of a3 + b1, a3 + b2, a4 + b2, and a1 + b3 in five-node combined networks and a3 + b1 and a3 + b4 in four-node combined networks. We find that the networks containing the hybrid motif a3 + b1 have the largest proportion; the mechanism for achieving the bi-function is given in Additional file 1. As three-node merging is full merging of the network, it is difficult to find very robust bi-functional networks with large Q1,2 values (no more than 5). Due to regulation conflicts, the F1 motifs (a1, a2, a3, and a4) and F2 motifs (b1, b2, b3, and b4) cannot be fully merged to form a bifunctional motif. Instead, the bi-functional networks of combination are derived from non-core motifs that can achieve both oscillation and adaptation, such as the structure skeletons c1 and c2 depicted in Fig. 3c.
By combining the nodes of functional F1 and F2 networks, the merged nodes function as a joint uniting the sub-functional networks. This results in competing effects between substrates and additional implicit interactions within the networks. This side effect of node merging could promote, hamper, or have little influence on the circuit’s function, depending on the strength of the competing strength. In Fig. 3d-f, we list several examples of bi-functional networks that contain implicit interactions as denoted by the dashed linkages. It is noted that the three-node bi-functional networks belong to the intersection of three-node networks in F1 and F2 pools, and that only a portion of networks in the intersection can achieve the bi-function. The node merging procedure does not discover new three-node networks that are bifunctional. For instance, all compatibly combined networks of oscillation motifs a1-a4 and adaptation motifs b1-b4 by full merging of the three nodes have a zero Q1,2 value.
Extending the bi-functional networks into tri-functional networks
Next, we extend the bi-functional networks to the target tri-functional networks. As illustrated in Fig. 4a, a bi-functional network can be processed in two ways to include an extra node (R3) for receiving the third signal. One approach is to adopt an existing regulatory node (i.e., a node other than R1, R2 and O) as the receptor for input of the third signal. The other strategy is to introduce a new node as R3 into a bifunctional network for receiving the third signal. As the third function requires sustained activation of the output node O, an activation from node R3 to the output node O is simply added. For the circuit in Fig. 4a, the node M1 or M2 can be used as the third receptor node R3; otherwise, an additional node R3 is added that activates the output node. Both strategies add no new competitive regulations to the networks.
For each candidate tri-functional network obtained from the above extension strategies, we sample 1000 random parameter sets for the newly introduced parameters while keeping the parameters in the original circuits unchanged. The number of parameter sets that can achieve the target tri-functional behavior is recorded as the Q(F3|F1,2)-value, which manifests the robustness for the triple signal-encoding function (see Additional file 1). The distribution of Q(F3| F1, 2) for all tri-functional networks obtained is shown in Fig. 4b. There are 3190 tri-functional networks with six nodes, 2092 tri-functional networks with five nodes, and 255 tri-functional networks with four nodes. The networks with Q(F3| F1, 2) > 500 are highlighted in the red frame in Fig. 4b, with 4 six-node networks, 134 five-node networks, and 3 four-node networks.
In Fig. 5c, we plot the structure diagram for all tri-functional networks with Q(F3| F1, 2) > 100 . A dot denotes a tri-functional network, and networks with Q(F3| F1, 2) > 500 are represented by large dots. All networks are distributed with respect to the number of links and nodes in the networks. We have linked any two networks that could be transformed by the gain or removal of one regulatory interaction. Most of the tri-functional networks we obtained are connected, and complex tri-functional networks can evolve from simpler core topologies. Several core networks are illustrated at the bottom in Fig. 5c, i.e., the networks labeled from N1 to N7. Network N4 is constructed by uniting motifs a3 and b1 and adding the R3 node. The network N3 is devised on the bi-functional networks with hybrid motifs of a1 + b3 via the first strategy. It has two possible nodes (R3 and M) used as a buffer node for F2 and can be reduced to the network N5 by removing node M. The robustness of network N5 is comparable with that of network N3. Thus, network N5 is the simplest robust tri-functional network in our calculation. In the following section, we perform an analysis of the dynamics of network N5.
Dynamics analysis for the simplest tri-functional network
For the simplest and functionally robust network N5 in Fig. 4c, we demonstrate its multiple signal-encoding dynamics in detail for illustrative purposes. The network dynamics in concentrations of active forms of enzymes R1, R2, R3 and O is described by the following ordinary differential equations,
In this description, the total concentrations of enzymes are normalized, and the variables vary from 0 to 1. In the normal condition, the receptor nodes R1, R2, R3 sense the surroundings with a base signal I1 = I2 = I3 = 0.1. With the parameter values in Fig. 5, the system rests at its stable steady state (with R1 = 0.1179, R2 = 0.3431, R3 = 0.0024, O = 0.1690). When node R1 is stimulated with I1 ≥ 0.2 and I2 = I3 = 0.1, the system enters into the regime of self-sustained oscillation (Fig. 5a). Oscillations are generated due to the R1-R3-O negative feedback loop within the circuit. The step signal is thus interpreted by the system into stable limit cycle oscillation in the output node. As depicted in Fig. 5d, the network dynamics undergoes a subcritical Hopf bifurcation at I1C = 0.25 when the strength of stimulus I1 is tuned as the control parameter.
For the second signal (I2 > 0.1, I1 = I3 = 0.1), the circuit’s output responds to the stimulus I2 transiently and returns to its pre-stimulated level. This transient activation results from the saturated regulations of R3 from both signal I3 and output node O, with K3 ≪ 1 − R3, K7 ≪ R3. In this case, Eq. 1c for R3 is simplified as,
At steady state, the above equation leads to the output which is proportional to the amplitude of I3,
The independence of OSS on I2 results in a perfect adaptation. Thus, the necessary condition for adaptation under the second stimulus requires that the node R3 functions as a buffer node within the negative feedback loop in the circuit. Figure 5b shows the responses of the output variable when node R2 is triggered by I2 of different strengths. The steady state of the network dynamics as a function of the signal amplitude I2 is depicted in Fig. 5e.
When node R3 of the circuit receives the third stimulus (I3 > 0.1, I1 = I2 = 0.1) of different amplitudes, the output node is always activated (see Fig. 5c), and the equilibrium state of the output increases with I3 strength (Fig. 5f). Therefore, the three signal-encoding processes are mutually compatible and completed in the simplest and robust topology.
Biochemical reactions in cellular signaling processes form complex molecular networks analogous to electrical circuits. There has been a long tradition of theoretical studies that looked into modules of small biochemical networks and signaling motifs that are viewed as simple building blocks of complex networks. Proteins functioning as the transfer and processing of information can be linked into motifs of reoccurring biochemical circuits that perform specific tasks such as to amplify signal, integrate and store information et al. [51, 52]. Motifs of common features that may help formalize the mapping from biochemical reactions to pathway block diagrams in signaling processes have been previously examined [53, 54]. In signaling and transcription regulation networks, the main classes of motifs that can carry out specific dynamic functions such as bistability, adaptation, oscillations, excitable pulse and et al. have been reviewed in references such as  and . Compared to previous studies that considered mainly mono-functional network motifs, we here report a systematic approach for extracting multi-functional motifs, i.e., depending on different external signals, the simple motifs can output dynamic responses of oscillation, transient activation, and sustained activation respectively.
The second novel point of this study lies in that competition between substrates for the same enzyme has been taken into account in our mathematical description of signaling circuits. Substrate competition that was usually omitted in previous theoretical studies can have pronounced effects. It has been reported that ultra-sensitivity can be generated through substrate competition between two sets of phosphorylation sites  or through negative cooperativity in which the binding of ligand to a multimeric receptor makes it more difficult for subsequent ligands to bind . Competitive binding can lead to high-dose hook effect (also known as prozone effect) in reactions where one protein acts as a linker between parts of a complex . In the interaction of Ca2+ and calcium binding proteins calmodulin (CaM), the competitive binding among CaM’s seven partners can tune the Ca2+/CaM binding frequency . In this study, we demonstrate that competitive effects could add implicit interactions within the networks that can promote or hamper the anticipated behaviors, and even create new dynamics unexpected from the explicit network topologies. For simple network structures without explicit functional motifs, the additional hidden regulations within the competitive nodes, combined with the rest of the linkages in the networks, could form implicit functional motifs that foster the networks to achieve target functions. For complex network topologies with dense linkages within the network, such as in the bi-functional networks obtained by two-node merging, the competitive effects could produce more complicated connections within the circuits that can drastically influence the network dynamics.
The multiple signal-encoding networks we consider involve at least four nodes. It is time-consuming to exhaust computationally all possible network topologies. We have adopted a compromise approach to efficiently construct multi-functional networks by enumerating modules of three nodes that can achieve each sub-function and then combing them by node merging. As the exploration is not exhaustive, the tri-functional networks we obtained may have occupied only a part of all functional topologies. As demonstrated, we have constructed tri-functional networks with four to six nodes via module combinations. The functions of oscillation and adaptation in the tri-functional networks could share a negative feedback loop, as long as the circuit could be driven to oscillate by the first signal and a node within as a buffer node for adaptation in the presence of the second signal. Alternatively, the function of adaptation could also be achieved via approaches such as the negative feedback loop with a buffer node (M2) in networks N1, N2, or N3. The third function of sustained activation, it requires at least one activation regulation directly or indirectly from node R3 to the output node O.
Compared to the functional networks we designed, the real signal-encoding network of Msn2 in Saccharomyces cerevisiae is much more complex [45,46,47,48,49,50]. A recent report that integrated experiments and modelling revealed that, the stimulus-dependent Msn2 dynamics is controlled by coupled feedback loops . At the core of this Msn2 pathway is the Ras–cAMP–PKA negative feedback loop which is crucial in shaping Msn2 dynamics. It is coupled to the signal-dependent mutual inhibition between protein kinase PKA and the yeast AMP-activated protein kinase Snf1. The negative feedback regulating PKA activity at the core of Msn2 pathway is abundant as sub-module in the functional networks that we identified here. As can be seen in Fig. 3 and Fig. 4, the negative feedback appears in forms of activator-inhibitor, delayed negative feedback, and repressilator. In the networks we found, coupled negative and positive feedback loops are not apparent. But when implicit interactions due to substrate competing are taken into account, such as those demonstrated in Fig. 2g, Fig. 2h, and in Fig. 3d, Fig. 3e, Fig. 3f, coupled negative and positive feedbacks are abundantly found as sub-modules in the networks identified here.
In p53 signaling pathway, the information of hundreds of cellular stress stimuli is encoded through the single node of p53 protein, and stimulus-specific p53 dynamics then triggers distinct cellular outcomes. While the real p53 signaling pathway is very complex, modelling studies  revealed that the flexibility of p53 activity is rooted in the core p53-Mdm2 negative feedback loop in the pathway: the transcription factor p53 activates the expression of Mdm2, and Mdm2 in turn inhibits p53 by either turning down its transcription or triggering its degradation. In the multifunctional networks that we find here, the negative feedback loop featuring the p53 pathway is also a fundamental sub-module. Negative feedback loops involving the output node are ubiquitous in bi- and tri- functional circuits we find here. In single cell study of p53 dynamics , γ-radiation and UV light cause DNA double strand breaks (DSBs) and single-stranded DNA (ssDNA) damage, respectively. It was uncovered that the mechanism of stimulus-specific response is relevant to two core coupled negative feedback loops. The upstream kinases ATM responding to DSBs and ATR responding to ssDNA both relay the damage signal to p53. This activates p53-Mdm2 and p53-Wip1 negative feedback loops. The core topology of two coupled negative feedback loops in p53 pathway in response to DNA damages is closely related to some of the networks that we identify here. As depicted in Fig. 3a, two coupled negative feedback loops are characteristic of the combined a3 + b2 and a4 + b2 networks, and apparently feature the combined networks in Fig. 3d and Fig. 3e. Although the tri-functional networks we constructed are simple and omit the real aspects of dynamic control, they could still be helpful for understanding the mechanism of real signal-encoding processes. Furthermore, the networks we obtained and the method we adopted could also be used in synthetic biology to construct multiple functional dynamic circuits for practical purposes.
The three signal-processing functions are compatible in a well-designed multi-functional network, in which the output dynamics is switched between oscillation, adaptation, and sustained activation by changing the input node for the signals, with no need to change the internal links in the networks. Competition between substrates for limited enzymes plays as a framework that renders implicit interactions between the nodes that are not explicitly linked. The substrate competition brings virtual links between nodes in the network and can create new implicit motifs that are seemingly nonfunctional.
We aim to search for networks capable of encoding the identity of distinct stimuli into the dynamics of the circuit’s output as a whole. The dynamics of oscillation, adaptation (transient activation) and sustained activation are the three target input-output responses considered. We assume that the functional network has three independent receptors each receiving a respective stimulus signal. Under normal conditions, all input signals are set as a low baseline, and the functional network remains in the resting state. When facing one of the three stimuli, the input signal is represented by a step function (Fig. 1a). The information about the stimulus is processed through the unknown components of the networks (Fig. 1b), and then the identity of the stimulus is encoded in the dynamics of the circuit’s output. The target dynamic outputs are illustrated in Fig. 1c as oscillations (F1), transient activation (F2), and sustained activation (F3). We intend to find the unknown networks in Fig. 1b that generate the target output when each receptor node receives its stimulus from an external signal.
In our calculations, the networks we consider are limited to enzymatic regulatory interactions. Each node in the network represents an enzymatic protein with a fixed total concentration (normalized to 1) and can be interconverted between its active and inactive forms by other enzymes that are active. Two types of links are considered. The arrow in the network denotes activation of the target enzyme, that is, the enzyme catalytically transforms its substrate from the inactive form to active form. The blunt-headed arrow indicates catalytic inactivation of the active substrate. In our model, we assume that nodes in the network can possibly be auto-activated or auto-inhibited. If a node in the network has no activation or deactivation from other nodes, it is assumed to be regulated by basally available enzymes in the background.
To take into account the competition effects among the enzymes that compete to bind the same target enzyme , we model the network dynamics with modified Michaelis-Menten rate equations. The dynamic variables represent the concentrations of enzymes (normalized to 1) in their active forms. For example, when enzyme R1 in the network converts both the enzymes M1 and O1 to active states M*1 and O*1 (Fig. 1d), the activation rate (RA) of M1 and O1 take the following forms (see additional file 1 for the full dynamical equations):
where R1, M1 and O1 are the concentrations of active enzymes R*1, M*1 and O*1, K1 and K2 denote the Michaelis-Menten constants, and V1 and V2 denote the regulation rate constants. Due to competition between enzymes M1 and O1 for binding enzyme R*1, an increase in the concentration of O*1 and therefore a decrease in inactive O1 (due to the constant total concentration) would increase the concentration of enzymes available for the inactive M1. This leads to an increased rate of M*1 production and vice versa. The situation is similar when enzyme R*1 deactivates both active enzymes M*1 and O*1. In these two cases, the competition effect results in implicit positive interactions between the two nodes as illustrated in Fig. 1d and Fig. 1e (note the dashed arrows in the circuits). Figure 1f illustrates the situation when enzyme R*1 activates M1 but deactivates O*1, resulting in implicit negative interactions of suppression. These competitive effects inducing implicit positive or negative interactions in enzymatic circuits have commonly been omitted in most theoretical or experimental models.
For a concrete demonstration of the competition effect between substrates, Fig. 1g and Fig. 1h demonstrate with a simple three-node network the distinct temporal behaviors generated with modified Michaelis-Menten equations and the normal Michaelis-Menten equations, respectively. With the same set of parameter values and same initial conditions, the description that takes into account competitive bindings depicts self-sustained oscillations (Fig. 1g) while a stationary steady state is generated when the normal Michaelis-Menten mechanism is adopted (Fig. 1h). The ordinary differential equations for the circuit in Fig. 1g and Fig.1 h are listed in additional file 1. One sees that the competition effect generates implicit mutual activation interactions as demonstrated with dashed links in Fig. 1g and sustains a completely different dynamical behavior. From the circuit that is added with two more edges, the oscillations not expected in the original topology are supported due to that a negative feedback loop forms with the aid of the implicit interactions coming from the competition effects. The competition effect could thus play significant roles in the topology of functional circuits and need to be considered in the design principles of multiple signal-encoding functional networks.
We adopt the method of module combination to construct the multi-functional networks. The scenario consists of three steps. First, we search separately for functional network topologies capable of oscillation and adaptation by enumerating all possible networks with three nodes. Typical and simple functional networks, i.e., those topologies that have less links and are relatively more robust, are then selected to form module pools for these two signal-encoding sub-functions. The combined network topologies are then constructed from the two sub-functional pools by merging nodes. Those that are capable of achieving both oscillation and adaptation in response to external stimuli are selected as bi-functional circuits. The final topologies for our target multiple signal-encoding functions are completed by adding to the bi-functional circuits an extra node that directly activates the output node upon receiving the signal, or alternatively by adopting an existing regulatory node for input of the third signal that directly or indirectly activates the output node. Ultimately, we obtain tri-functional networks having four to six nodes.
The modified Michaelis-Menten equations for network dynamics are integrated numerically. The code is written in computer C language with Visual Studio 2008. The ODE solver rkf45 of Runge-Kutta algorithm is chosen. The scripts are available upon request. Rate constant values ranging from 0.1 to 10 are considered. The Michaelis-Menten constant changes in the range 10−3~10. In order to determine computationally whether oscillations emerge out or not, the variance of network output variable in a time window is monitored after the transient behavior dies out. An apparently not zero variance indicates oscillations arise. The oscillations with both amplitude and oscillation period larger than 0.1 are considered in this study. To numerically determine adaptation, the transient peak value Opeak and asymptotic end value Oend of the output variable are monitored and calculate their differences from the initial value Oini. The numerical criteria for adaptation are that Opeak − Oini > 0.2 and ∣Oend − Oini ∣ < 0.1 are both satisfied. Sustained activation is determined by the criterion Oend − Oini > 0.2.
In order to analyze topological properties of three-node networks with oscillations and adaptation, we perform clustering analyses for the functional networks as in . In a three-node network, there are totally nine possible links because there are at most three links from each node to other nodes and to itself. Values of 1, − 1, and 0 are assigned for links of activation, inhibition, and no regulation, respectively. A network topology can be therefore described with a digital number of length nine. The topological difference between two networks having three nodes can then be measured by calculating the hamming pair-wise distance from the digital numbers. The clustering analyses of the functional networks is performed by using clustergram in Matlab. The clustering for bi-functional networks is similarly analyzed.
The original c code files used in our computation for enumeration of sub-functions (Additional file 2), and merging of functional modules (Additional files 3, 4, 5) as well as a model in SBML format (Additional file 6) for the typical network (N5 in Fig. 4c) can be found in Additional files.
Jimenez A, Cotterell J, Munteanu A, Sharpe J. A spectrum of modularity in multi-functional gene circuits. Mol Syst Biol. 2017;13(4):925.
Lim WA, Lee CM, Tang C. Design principles of regulatory networks: searching for the molecular algorithms of the cell. Mol Cell. 2013;49(2):202–12.
Purvis JE, Lahav G. Encoding and decoding cellular information through signaling dynamics. Cell. 2013;152:954–6.
Hao N, Budnik BA, Gunawardena J, O'Shea EK. Tunable signal processing through modular control of transcription factor translocation. Science. 2013;339(6118):460–4.
Levine JH, Lin Y, Elowitz MB. Functional roles of pulsing in genetic circuits. Science. 2013;342(6163):1193–200.
Hao N, O'Shea EK. Signal-dependent dynamics of transcription factortranslocation controls gene expression. Nat Struct Mol Biol. 2012;19:31–40.
Locke JCW, Young JW, Fontes M, Jiménez MJH, Elowitz MB. Stochastic pulse regulation in bacterial stress response. Science. 2011;334(6054):366.
Sadeh A, Movshovich N, Volokh M, Gheber L, Aharoni A. Fine-tuning of the Msn2/4–mediated yeast stress responses as revealed by systematic deletion of Msn2/4 partners. Mol Biol Cell. 2011;22(17):3127–38.
Mitchell A, Wei P, Wendell AL. Oscillatory stress stimulation uncovers an Achilles' heel of the yeast MAPK signaling network. Science. 2015;350:1379–83.
AkhavanAghdam Z, Sinha J, Tabbaa OP, Hao N. Dynamic control of gene regulatory logic by seemingly redundant transcription factors. eLife. 2016;5:e18458.
Li Y, Roberts J, AkhavanAghdam Z, Hao N. Mitogen-activated protein kinase (MAPK) dynamics determine cell fate in the yeast mating response. J Biol Chem. 2017;292:20354.
Benzinger D, Khammash M. Pulsatile inputs achieve tunable attenuation of gene expression variability and graded multi-gene regulation. Nat Commun. 2018;9:3521.
Nandagopal N, Santat LA, LeBon L, Sprinzak D, Bronner ME, Elowitz MB. Dynamic ligand discrimination in the notch signaling pathway. Cell. 2018;172:869–80.
Selimkhanov J, Taylor B, Yao J, Pilko A, Albeck J, Hoffmann A, et al. Accurate information transmission through dynamic biochemical signaling networks. Science. 2014;346(6215):1370.
O'Shaughnessy EC, Palani S, Collins JJ, Sarkar CA. Tunable signal processing in synthetic MAP kinase cascades. Cell. 2011;144(1):119–31.
Hansen AS, O'Shea EK. Promoter decoding of transcription factor dynamics involves a trade-off between noise and control of gene expression. Mol Syst Biol. 2013;9:704.
Batchelor E, Loewer A, Mock C, Lahav G. Stimulus-dependent dynamics of p53 in single cells. Mol Syst Biol. 2011;7(1):488.
Purvis JE, Karhohs KW, Mock C, Batchelor E, Loewer A, Lahav G. p53 dynamics control cell fate. Science. 2012;336(6087):1440–4.
Garmendia-Torres C, Goldbeter A, Jacquet M. Nucleocytoplasmic oscillations of the yeast transcription factor Msn2: evidence for periodic PKA activation. Curr Biol. 2007;17(12):1044–9.
De Wever V, Reiter W, Ballarini A, Ammerer G, Brocard C. A dual role for PP1 in shaping the Msn2-dependent transcriptional response to glucose starvation. EMBO J. 2005;24(23):4115–23.
Hohmann S. Osmotic stress signaling and Osmoadaptation in yeasts. Microbiol Mol Biol Rev. 2002;66(2):300–72.
Garcia R, Rodriguez-Pena JM, Bermejo C, Nombela C, Arroyo J. The high osmotic response and cell wall integrity pathways cooperate to regulate transcriptional responses to zymolyase-induced cell wall stress in Saccharomyces cerevisiae. J Biol Chem. 2009;284(16):10901–11.
Mettetal JT, Muzzey D, Gomez-Uribe C, van Oudenaarden A. The frequency dependence of osmo-adaptation in Saccharomyces cerevisiae. Science. 2008;319(5862):482–4.
Gat-Viks I, Shamir R. Refinement and expansion of signaling pathways: the osmotic response network in yeast. Genome Res. 2007;17(3):358–67.
Gennemark P, Nordlander B, Hohmann S, Wedelin D. A simple mathematical model of adaptation to high osmolarity in yeast. In Silico Biol. 2006;6(3):193–214.
Charizanis C, Juhnke H, Krems B, Entian KD. The oxidative stress response mediated via Pos9/Skn7 is negatively regulated by the Ras/PKA pathway in Saccharomyces cerevisiae. Mol Gen Genet. 1999;261(4–5):740–52.
Hasan R, Leroy C, Isnard AD, Labarre J, Boy-Marcotte E, Toledano MB. The control of the yeast H2O2 response by the Msn2/4 transcription factors. Mol Microbiol. 2002;45(1):233–41.
Tay S, Hughey JJ, Lee TK, Lipniacki T, Quake SR, Covert MW. Single-cell NF-κB dynamics reveal digital activation and analog information processing in cells. Nature. 2010;466(7303):267.
Sasagawa S, Y-i O, Fujita K, Kuroda S. Prediction and validation of the distinct dynamics of transient and sustained ERK activation. Nat Cell Biol. 2005;7(4):365–73.
Huang S, Ingber DE. Shape-dependent control of cell growth, differentiation, and apoptosis: switching between attractors in cell regulatory networks. Exp Cell Res. 2000;261(1):91–103.
Wagner A. Circuit topology and the evolution of robustness in two-gene circadian oscillators. Proc Natl Acad Sci U S A. 2005;102(33):11775–80.
Castillo-Hair SM, Villota ER, Coronado AM. Design principles for robust oscillatory behavior. Syst Synth Biol. 2015;9(3):125–33.
Kobayashi Y, Shibata T, Kuramoto Y, Mikhailov AS. Evolutionary design of oscillatory genetic networks. Euro Phys J B. 2010;76(1):167–78.
Shah NA, Sarkar CA. Robust network topologies for generating switch-like cellular responses. PLoS Comput Biol. 2011;7(6):e1002085.
Zeng W, Du P, Lou Q, Wu L, Zhang HM, Lou C, et al. Rational Design of an Ultrasensitive Quorum-Sensing Switch. ACS Synth Biol. 2017;6(8):1445–52.
Ma W, Trusina A, El-Samad H, Lim WA, Tang C. Defining network topologies that can achieve biochemical adaptation. Cell. 2009;138(4):760–73.
Shi W, Ma W, Xiong L, Zhang M, Tang C. Adaptation with transcriptional regulation. Sci Rep. 2017;7:42648.
Paul F, Eric DS. A case study of evolutionary computation of biochemical adaptation. Phys Biol. 2008;5(2):026009.
Yan L, Ouyang Q, Wang H. Dose-response aligned circuits in signaling systems. PLoS One. 2012;7(4):e34727.
Xi JY, Ouyang Q. Using sub-network combinations to scale up an enumeration method for determining the network structures of biological functions. PLoS One. 2016;11(12):e0168214.
Bintu L, Buchler NE: Garcia HG, Gerland U, Hwa T, Kondev J, et al.: Transcriptional regulation by the numbers: models. Curr Opin Genet Dev 2005;15(2):116–124.
Rondelez Y. Competition for catalytic resources alters biological network dynamics. Phys Rev Lett. 2012;108(1):018102.
Iman RL, Davenport JM, Zeigler DK: Latin hypercube sampling (program user's guide).[LHC, in FORTRAN]. Sandia Labs., Albuquerque, NM (USA); 1980.
Cotterell J, Sharpe J. An atlas of gene regulatory networks reveals multiple three-gene mechanisms for interpreting morphogen gradients. Mol Syst Biol. 2010;6:425.
Gutin J, Sadeh A, Rahat A, Aharoni A, Friedman N. Condition-specific genetic interaction maps reveal crosstalk between the cAMP/PKA and the HOG MAPK pathways in the activation of the general stress response. Mol Syst Biol. 2015;11(10):829.
Chasman D, Ho YH, Berry DB, Nemec CM, MacGilvray ME, Hose J, et al. Pathway connectivity and signaling coordination in the yeast stress-activated signaling network. Mol Syst Biol. 2014;10(11):759.
Jacquet M, Renault G, Lallet S, De Mey J, Goldbeter A. Oscillatory nucleocytoplasmic shuttling of the general stress response transcriptional activators Msn2 and Msn4 in Saccharomyces cerevisiae. J Cell Biol. 2003;161(3):497–505.
Natalia Petrenkoa RzVC, Megan N. McCleanc, Alexandre V. Morozovb,d, James R. Broacha: Noise and interlocking signaling pathways promote distinct transcription factor dynamics in response to different stresses. Mol Biol Cell 2013;24:2045–2057.
Morano KA, Grant CM, Moye-Rowley WS. The response to heat shock and oxidative stress in Saccharomyces cerevisiae. Genetics. 2012;190(4):1157–95.
Boisnard S, Lagniel G, Garmendia-Torres C, Molin M, Boy-Marcotte E, Jacquet M, et al. H2O2 activates the nuclear localization of Msn2 and Maf1 through thioredoxins in Saccharomyces cerevisiae. Eukaryot Cell. 2009;8(9):1429–38.
Bray D. Protein molecules as computational elements in living cells. Nature. 1995;376:307–12.
Bhalla US, Iyengar R. Emergent properties of networks of biological signaling pathways. Science. 1999;283:381–7.
Bhalla US. The chemical organization of signaling interactions. Bioinformatics. 2002;18:855–63.
Tyson JJ, Chen KC, Novak B. Sniffers, buzzers, toggles and blinkers: dynamics of regulatory and signaling pathways in the cell. Curr Opin Cell Biol. 2003;15(2):221–31.
Shoval O, Alon U. SnapShot: network motifs. Cell. 2010;143:326–e1.
Kim SY, Ferrell JJE. Substrate competition as a source of ultrasensitivity in the inactivation of Wee1. Cell. 2007;128:1133–45.
Ha SH, Ferrell JJE. Thresholds and ultrasensitivity from negative cooperativity. Science. 2016;352:990–3.
Roy RD, Rosenmund C, Stefan MI. Cooperative binding mitigates the high-dose hook effect. BMC Syst Biol. 2017;11:74.
Romano DR, Pharris MC, Patel NM, Kinzer-Ursem TL. Competitive tuning: Competition’s role in setting the frequency- dependence of Ca2+-dependent proteins. PLoS Comput Biol. 2017;13(11):e1005820.
Jiang Y, AkhavanAghdam Z, Tsimring ALS, Hao N. Coupled feedback loops control the stimulus-dependent dynamics of the yeast transcription factor Msn2. J Biol Chem. 2017;292:12366–72.
Hunziker A, Jensen MH, Krishna S. Stress-specific response of the p53-Mdm2 feedback loop. BMC Syst Biol. 2010;4:94.
The authors thank the financial support of the Ministry of Science and Technology of China (2015CB910300 to HW) and National Natural Science Foundation of China (11434001 to QO).
HW was supported by the Ministry of Science and Technology of China (2015CB910300). QO was supported by the National Natural Science Foundation of China (11434001).
Availability of data and materials
The codes and SBML file for generating the data are available in Additional files.
Consent to publish
Ethics approval and consent to participate
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
The appendix for the equations and parameters used in Fig.1(d-f), and module selection for the function of sustained activation, as well as the analysis of bi-functional mechanism for the hybrid motif of a3 + b1. (DOCX 165 kb)
The original codes used in our computation for enumeration of sub-functions. (ZIP 2318 kb)
The original codes used in our computation for one-node merging of functional modules. (ZIP 1608 kb)
The original codes used for two-node merging of functional modules. (ZIP 1018 kb)
The original codes used for three-node merging of functional modules. (ZIP 1162 kb)