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 [43]. 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 [31], 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 [36], 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 [36]. 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:
-
(i).
Output O nodes have to be merged with each other as we assume a unique output node.
-
(ii).
M nodes may be merged with each other (unless there are contradictory linkages).
-
(iii).
R nodes may be merged with M nodes from the other pathway (unless there are contradictory linkages).
-
(iv).
R nodes for receiving signals cannot be merged with each other.
-
(v).
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 [44]. 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,
$$ \left\{\begin{array}{c}\begin{array}{c}\frac{d{R}_1}{dt}={V}_1{I}_1\frac{1-{R}_1}{K_1+1-{R}_1}-{V}_5{R}_3\frac{R_1}{K_5+{R}_1}\\ {}\frac{d{R}_2}{dt}={V}_2{I}_2\frac{1-{R}_2}{K_2+1-{R}_2}-{V}_6B\frac{R_2}{K_6+{R}_2}\end{array}\\ {}\begin{array}{c}\frac{d{R}_3}{dt}={V}_3{I}_3\frac{1-{R}_3}{K_3+1-{R}_3}-{V}_7O\frac{R_3}{K_7+{R}_3}\\ {}\frac{dO}{dt}={V}_4{R}_2\frac{1-O}{K_4+1-O}-{V}_8{R}_1\frac{O}{K_8+O}\end{array}\end{array}\right. $$
(1)
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,
$$ \frac{d{R}_3}{dt}\approx {V}_3{I}_3-{V}_7O $$
(2)
At steady state, the above equation leads to the output which is proportional to the amplitude of I3,
$$ {O}_{SS}=\frac{V_3{I}_3}{V_7} $$
(3)
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.