Module selection for each signalprocessing function
We use the method of threenodenetwork enumeration to explore the functional networks of oscillation (F_{1}) and adaptation (F_{2}), respectively. In a two or threenode 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 subfunctions F_{1} and F_{2} is measured by the number of parameter sets, i.e., Q_{1} and Q_{2}values, that generate the target dynamic functions. For functions F_{1} and F_{2}, we select the network topologies with relatively high Qvalues and fewer edges as functional F_{1} and F_{2} module pools.
The oscillation (F_{1}) pool contains 73 networks with 3 to 6 edges with Q_{1} > 50 and 8 simplest networks with 3 edges with 10 < Q_{1} < 50 (as marked in Fig. 2a). Clustering analysis of these networks (Fig. 2b) shows that the basic oscillation motifs for achieving function F_{1} are negative feedback loops between the three nodes. Motifs a_{1} and a_{3} are simply two forms of the repressilator, and motifs a_{2} and a_{4} are the delayed feedback loop. All 73 networks contain explicitly one of the two negative feedback loops which are wellknown motifs for oscillations. This is true even for networks with edge numbers greater than 6. Additionally, the F_{1} 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 selfsustained 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 a_{3} and a_{4} motifs. Thus, the modular pool of F_{1} consists of 81 networks as marked in Fig. 2a.
Similarly, as marked in Fig. 2d, the F_{2} functional pool also consists of 81 selected networks, i.e., 76 networks with 3 to 6 edges and Q_{2} > 55, and 5 networks with 3 edges and 10 < Q_{2} < 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 b_{1}, b_{3}), a delayed feedback loop (motif b_{2}), and an incoherent feedforward loop with a proportional node (motif b_{4}). The motifs for function F_{2} are similar to Ma’s results [31], but different in that motifs b_{3} and b_{4} (for which Q_{2} < 55) are less robust than motifs b_{1} and b_{2} (for which Q_{2} > 55). This difference is attributable to the inclusion of competitive regulations in both motifs b_{3} and b_{4}: node M_{3} and node O_{3} compete for enzyme R_{3} to convert from the inactive to active state. Consequently, there exit implicit positive interactions between nodes M_{3} and O_{3}, 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 Q_{2} value as marked in Fig. 2d. These simple functional circuits all have implicit interactions, of which two are motifs b_{3} and b_{4}. 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 b_{1} and b_{3}.
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 (F_{3}) are not searched by enumeration. We prefer to adopt the simplest structure, i.e., a link of activation from the input node R_{3} to the output node O_{3}. Thus, we have constructed the module pools for the subfunctions of oscillation (F_{1}, with 81 topologies), transient activation (F_{2}, with 81 topologies), and sustained activation (F_{3}, with 1 topology), respectively.
Module combination for bifunctional networks of F_{1} and F_{2}
By combining subfunction networks in the F_{1} and F_{2} module pools, we next construct F_{1}F_{2} bifunctional networks that encode the external signal received by nodes R_{1} and R_{2} respectively into oscillations and adaptation in the output node O. The network topologies from the two subfunctional pools are united by merging their nodes. As demonstrated in Fig. 3, we adopt the following nodemerging 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 bifunctional network, there are three possible ways to combine the F_{1} and F_{2} circuits as illustrated in the upper left of Fig. 3ac. One approach is to simply merge the output nodes O_{1} and O_{2}. Another is twonode merging of M_{1}M_{2}, R_{1}M_{2}, and R_{2}M_{2} in addition to the O_{1}O_{2} merging. The third approach is the threenode merging of O_{1}O_{2}, R_{1}M_{2}, and R_{2}M_{1}.
We enumerate all possible networks combined from the two module pools as candidate topologies for bifunctional 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 Q_{1,2}value to estimate the robustness of the bifunction. In the bifunctional networks with Q_{1,2} > 0, we obtain 3215 fivenode, 1337 fournode, and 41 threenode networks by one, two, and threenode merging, respectively. The Q_{1,2} distributions of these three classes of bifunctional combined networks are illustrated in the lower left of Fig. 3ac. The networks with relatively high Q_{1,2} are highlighted with red dashed boxes in the histograms of Fig. 3. These robust networks obtained by one, two and threenode merging are then clustered (middle panels in Fig. 3ac). Structure skeletons from the clustering analyses are listed alongside (right panels in Fig. 3ac). The robust bifunctional networks obtained via onenode and twonode merging are composed of hybrid motifs from the F_{1} and F_{2} module pools, such as hybrids of a_{3} + b_{1}, a_{3} + b_{2}, a_{4} + b_{2}, and a_{1} + b_{3} in fivenode combined networks and a_{3} + b_{1} and a_{3} + b_{4} in fournode combined networks. We find that the networks containing the hybrid motif a_{3} + b_{1} have the largest proportion; the mechanism for achieving the bifunction is given in Additional file 1. As threenode merging is full merging of the network, it is difficult to find very robust bifunctional networks with large Q_{1,2} values (no more than 5). Due to regulation conflicts, the F_{1} motifs (a_{1}, a_{2}, a_{3}, and a_{4}) and F_{2} motifs (b_{1}, b_{2}, b_{3}, and b_{4}) cannot be fully merged to form a bifunctional motif. Instead, the bifunctional networks of combination are derived from noncore motifs that can achieve both oscillation and adaptation, such as the structure skeletons c_{1} and c_{2} depicted in Fig. 3c.
By combining the nodes of functional F_{1} and F_{2} networks, the merged nodes function as a joint uniting the subfunctional 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. 3df, we list several examples of bifunctional networks that contain implicit interactions as denoted by the dashed linkages. It is noted that the threenode bifunctional networks belong to the intersection of threenode networks in F1 and F2 pools, and that only a portion of networks in the intersection can achieve the bifunction. The node merging procedure does not discover new threenode networks that are bifunctional. For instance, all compatibly combined networks of oscillation motifs a_{1}a_{4} and adaptation motifs b_{1}b_{4} by full merging of the three nodes have a zero Q_{1,2} value.
Extending the bifunctional networks into trifunctional networks
Next, we extend the bifunctional networks to the target trifunctional networks. As illustrated in Fig. 4a, a bifunctional network can be processed in two ways to include an extra node (R_{3}) for receiving the third signal. One approach is to adopt an existing regulatory node (i.e., a node other than R_{1}, R_{2} and O) as the receptor for input of the third signal. The other strategy is to introduce a new node as R_{3} 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 R_{3} to the output node O is simply added. For the circuit in Fig. 4a, the node M_{1} or M_{2} can be used as the third receptor node R_{3}; otherwise, an additional node R_{3} is added that activates the output node. Both strategies add no new competitive regulations to the networks.
For each candidate trifunctional 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 trifunctional behavior is recorded as the Q(F_{3}F_{1,2})value, which manifests the robustness for the triple signalencoding function (see Additional file 1). The distribution of Q(F_{3} F_{1, 2}) for all trifunctional networks obtained is shown in Fig. 4b. There are 3190 trifunctional networks with six nodes, 2092 trifunctional networks with five nodes, and 255 trifunctional networks with four nodes. The networks with Q(F_{3} F_{1, 2}) > 500 are highlighted in the red frame in Fig. 4b, with 4 sixnode networks, 134 fivenode networks, and 3 fournode networks.
In Fig. 5c, we plot the structure diagram for all trifunctional networks with Q(F_{3} F_{1, 2}) > 100 [44]. A dot denotes a trifunctional network, and networks with Q(F_{3} F_{1, 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 trifunctional networks we obtained are connected, and complex trifunctional networks can evolve from simpler core topologies. Several core networks are illustrated at the bottom in Fig. 5c, i.e., the networks labeled from N_{1} to N_{7}. Network N_{4} is constructed by uniting motifs a_{3} and b_{1} and adding the R_{3} node. The network N_{3} is devised on the bifunctional networks with hybrid motifs of a_{1} + b_{3} via the first strategy. It has two possible nodes (R_{3} and M) used as a buffer node for F_{2} and can be reduced to the network N_{5} by removing node M. The robustness of network N_{5} is comparable with that of network N_{3}. Thus, network N_{5} is the simplest robust trifunctional network in our calculation. In the following section, we perform an analysis of the dynamics of network N_{5}.
Dynamics analysis for the simplest trifunctional network
For the simplest and functionally robust network N_{5} in Fig. 4c, we demonstrate its multiple signalencoding dynamics in detail for illustrative purposes. The network dynamics in concentrations of active forms of enzymes R_{1}, R_{2}, R_{3} 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{1O}{K_4+1O}{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 R_{1}, R_{2}, R_{3} sense the surroundings with a base signal I_{1} = I_{2} = I_{3} = 0.1. With the parameter values in Fig. 5, the system rests at its stable steady state (with R_{1} = 0.1179, R_{2} = 0.3431, R_{3} = 0.0024, O = 0.1690). When node R_{1} is stimulated with I_{1} ≥ 0.2 and I_{2} = I_{3} = 0.1, the system enters into the regime of selfsustained oscillation (Fig. 5a). Oscillations are generated due to the R_{1}R_{3}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 I_{1C} = 0.25 when the strength of stimulus I_{1} is tuned as the control parameter.
For the second signal (I_{2} > 0.1, I_{1} = I_{3} = 0.1), the circuit’s output responds to the stimulus I_{2} transiently and returns to its prestimulated level. This transient activation results from the saturated regulations of R_{3} from both signal I_{3} and output node O, with K_{3} ≪ 1 − R_{3}, K_{7} ≪ R_{3}. In this case, Eq. 1c for R_{3} 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 I_{3},
$$ {O}_{SS}=\frac{V_3{I}_3}{V_7} $$
(3)
The independence of O_{SS} on I_{2} results in a perfect adaptation. Thus, the necessary condition for adaptation under the second stimulus requires that the node R_{3} functions as a buffer node within the negative feedback loop in the circuit. Figure 5b shows the responses of the output variable when node R_{2} is triggered by I_{2} of different strengths. The steady state of the network dynamics as a function of the signal amplitude I_{2} is depicted in Fig. 5e.
When node R_{3} of the circuit receives the third stimulus (I_{3} > 0.1, I_{1} = I_{2} = 0.1) of different amplitudes, the output node is always activated (see Fig. 5c), and the equilibrium state of the output increases with I_{3} strength (Fig. 5f). Therefore, the three signalencoding processes are mutually compatible and completed in the simplest and robust topology.