Entangled gene regulatory networks with cooperative expression endow robust adaptive responses to unforeseen environmental changes

Living organisms must respond to environmental changes. Generally, accurate and rapid responses are provided by simple, unidirectional networks that connect inputs with outputs. Besides accuracy and speed, biological responses should also be robust to environmental or intracellular noise and mutations. Furthermore, cells must also respond to unforeseen environmental changes that have not previously been experienced, to avoid extinction prior to the evolutionary rewiring of their networks, which takes numerous generations. We have investigated gene regulatory networks that mutually activate or inhibit, and have demonstrated that complex entangled networks can make appropriate input-output relationships that satisfy the robust and adaptive responses required for unforeseen challenges. Such entangled networks function for sloppy and unreliable responses with low Hill coefficient reactions for the expression of each gene. To compensate for such sloppiness, several detours in the regulatory network exist. By taking advantage of the averaging over such detours, the network shows a higher robustness to environmental and intracellular noise as well as to mutations in the network, when compared to simple unidirectional circuits. Furthermore, the appropriate response to unforeseen challenges, allowing for functional outputs, is achieved as many genes exhibit similar dynamic expression responses, irrespective of inputs, as confirmed by applying dynamic time warping and dynamic mode decomposition. As complex entangled networks are common in gene regulatory networks and global gene expression responses are observed in microbial experiments, the present results provide a novel design principle for cellular networks.

Recent experimental advances have demonstrated that cells often have appropriate, robust responses to environmental changes, including those that have not previously been experienced. It is known that accurate and rapid responses can be achieved by simple unidirectional networks that connect straightforwardly input and outputs. However, such responses were not robust to perturbations. Here we have made numerical evolution of gene regulatory networks with mutual activation and inhibitions, and uncovered that complex entangled networks including many feedforward and feedback paths can make robust input-output responses, when each gene expression is not accurate. Remarkably, they make appropriate responses even to unforeseen environmental changes, as are supported by global, correlated responses across genes that are similar for all input signals. The results explain why cells adopt complex gene regulatory networks and exhibit global expression changes, even though they may not be advantageous in terms of their energy cost or response speed. The present results are consistent with the recent experiments on microbial gene expression changes and network analyses. This investigation provides insights into how cells survive fluctuating and unforeseen unpredictable environmental changes, and gives a universal conceptual framework to go beyond the standard picture based on a combination of network motifs.

Introduction
Living organisms, generally respond appropriately to environmental changes by adapting their internal states to the new conditions. In cells, chemical compositions, in particular gene expression levels, are adjusted to adapt to novel environments. Information on the environmental change is transmitted via a signal transduction network to regulate gene expression patterns.
For a cell to make an appropriate response rapidly and accurately, it will be relevant to transfer the input signal to the output expression pattern unidirectionally through a sequence of chemical responses. Indeed, such unidirectional circuits are ubiquitously observed as direct connections or feed-forward networks from the input in the signal transduction networks. How these unidirectional networks achieve appropriate input-output relationships has previously been extensively investigated [1][2][3].
Real cellular networks, however, are often much more complex. They include many components and are entangled. Sometimes such unidirectional circuits may be extracted as a part of the complex network, but it is not clear whether they have the same function as has been determined in isolation [4,5]. They generally do not function independently, as the networks are entangled. Why cells have adopted such complex networks, however, is not completely understood.
There are further requirements for cellular responses beyond being quick and accurate. For instance, they should be robust to perturbations [6][7][8][9][10], as the environmental conditions surrounding living organisms normally fluctuate [11][12][13]. Moreover, each process controlling the cellular concentrations of chemicals is generally stochastic, as the number of molecules involved is sometimes not very large [14][15][16]. In addition to the robustness to fluctuations and noise, there is another postulate, an 'adaptive response to an unforeseen challenge (AUC)', as put forward by Braun [17]: Cells must manage to respond to and survive under unforeseen environmental changes that are not expected by the prescribed input-output responses. In particular, even when cells experience a novel environment for which appropriate output responses are not yet prepared, they have to survive to some degree to avoid extinction, as it takes many generations for the appropriate evolutionary changes to occur. Responses to unforeseen challenges are thus required without rewiring the gene regulatory network (GRN) via genetic evolution [18][19][20].
Direct unidirectional networks optimized for given input-output relationships are not expected to be robust. Further, it will be difficult to achieve AUCs with them. If this is the case, and if a complex network can have advantages over unidirectional networks, with regard to robustness and AUCs, this may explain why biological networks are so complex and entangled. However, it is not yet evident if a complex network with many degrees of freedom acquires such robustness and AUCs, as they do not always have functional robustness. It is thus necessary to understand which types of complex networks can have a higher capacity for robustness and AUCs than simple unidirectional networks and to uncover design principles for such networks.
Previously we have evolved GRNs to achieve appropriate input-output relationships. In addition to simple direct or feed-forward networks, we have uncovered another type of entangled network consisting of many components, that can generate an appropriate input-output relationship as a cooperative response of many genes [21]. In particular, this cooperative response emerges when each gene expression is sloppy, that is, with a low Hill coefficient as observed in a real GRN [22][23][24][25]. Here, we explore whether this class of networks achieves higher capacity in noise robustness and AUCs than traditional direct or feed-forward networks.

Results
We adopted a simplified GRN model [8,21,[26][27][28][29][30] and investigated the networks evolved so that, for each given input, one prescribed output gene (target) is expressed, by defining the fitness of the expression patterns of the given output genes (Eq (2)). Each input stimulated one of the input genes and was transmitted to the output genes through the middle-layer genes (ML-genes). As described in [21], three types of networks evolved depending on the sensitivity (corresponding to the Hill coefficient) of the gene expression dynamics: direct connections (Direct type) and feed-forward networks having side paths (FF-network type) that showed rapid and accurate responses evolved when each gene's response was sensitive (large β), whereas complex entangled networks with cooperative responses of many genes (Cooperative type) evolved when each response was sloppy (small β). The three types showed different characteristics in their network structures, and these were most discernible when their core structures were compared. The core structure was obtained by removing paths successively, one-by-one, as long as the corresponding output genes response for each input was preserved (e.g., maintaining fitness > 0.8) (Fig 1). The Direct and FF-network types adopted unidirectional circuits from the input to the output genes. The Cooperative network type, however, involved many genes and its network was entangled. In this case, with the input, a few genes were locally excited, whereas global inhibition by many other genes followed. The characteristic structure of local excitation and global inhibition (LEGI [31]) was common in this Cooperative type. Here, it is of note that the evolutionary simulations were carried out to preserve the total path number of the networks, and the redundant paths out of the core structure remained for all types.

Noise robustness
Robustness to noise was investigated first. A white Gaussian noise term with amplitude ξ was added to the time evolution of each expression level (Eq (1)) for the ML-genes (Fig 2 (a)) or stimulated input genes (Fig 2  (b)), and the corresponding Langevin equations were simulated. When the noise term was added to the MLgenes, the Direct type was the most fragile and the fitness was drastically reduced by the noise (Fig 2 (a)). The FF-network and the Cooperative types showed noise robustness, as their signal transduction routes were multitracked and noise could be canceled out along the pathways. Furthermore, the Cooperative type was more robust to smaller amounts of noise, with ξ < 0.1. This was because the entangled network enabled the averaging out of noise, not only on each gene at the endpoint of multitrack paths but also on every gene in the network. It is of note that the noise canceling effect was independent of how many genes affected each output gene. An output gene had interaction paths with 6.1 genes in the Direct type, 13.0 genes in the FF-network type, and 11.0 genes in the Cooperative type, on average. In addition, the total number of paths in a network was set to be the same for all types. Hence, there was no correlation between the number of paths of interacting genes and noise robustness.
We also studied robustness against noise to an input gene (Fig 2 (b)). Again, the Direct type was largely influenced by the input noise as the noise effect was transmitted directly via the straight paths. The FF-network type was less robust than the Direct type. This is because its feed-forward network amplifies the input noise as well as the external signal. The Cooperative type had higher robustness to the input noise. Again, the noise could be reduced via the detour paths in the pathways.
Here, the path distance between input and output genes (i.e., the number of intermediate ML-genes between them) was irrelevant to the noise robustness. The shortest distance was similar for the three types, and the input and corresponding output genes were connected by one or two ML-genes. The noise was reduced not through a long-connected pathway but through multitrack, parallel paths.

Mutational robustness
Robustness against the mutational changes in the GRN were then considered. Here, some paths, especially those that constituted a core structure, were inevitable for responses, but some other paths were not necessary to achieve high levels of fitness. The robustness against the removal of the latter was trivial. Thus, mutational robustness should be defined as a fitness change when indispensable paths were removed. Consequently, we studied mutational robustness against the number of removed paths that constituted a core structure (Fig 2 (c)).
The Direct type had no mutational robustness, as the response was destroyed when any path that constituted a core structure was removed. The other two network types were more robust against mutational changes as their core structures contained detour paths. Even when a route connecting an input-output pair was disconnected, the external signal could still be transmitted by way of the detour paths. However, the FF-network and the Cooperative types were different in that the core structure of the former was separated according to each input-output pair, but not in the latter. Hence, their mechanisms of mutational robustness differed. In the FF-network type, when a path forming a core structure was removed, the response of a corresponding output gene dropped off, but the responses of the other output genes were hardly influenced. In contrast, in the Cooperative type with the LEGI structure, the responses of all output genes were weakened when any path in a core structure were removed. However, the decline was smaller as the Cooperative type included a larger number of alternative detour paths.
In the evolutionary process so far, we have fixed each parameter (response sensitivity corresponding to the Hill coefficient β and a constant threshold for expression y T ) of each gene to a common value and optimized the network structure for the given parameter values. We also studied the robustness to the changes in these parameter values (Fig 3). The Direct and FF-network types maintained a response only within a parameter range for each type, and the response often disappeared outside of their parameter range. On the other hand, the Cooperative type could maintain the response beyond the parameter range corresponding to its phase. This is because the LEGI structure includes the core structures for both the Direct and FF-network types within.

Tolerance to erroneous input
Thus far, the external input S ext was maintained at a steady value for t ≥ 0 in our model. Responses to a pulse input signal that was applied only for a brief time (0 ≤ t ≤ τ ) were then investigated. These signals represent an erroneous input due to environmental fluctuations. It is of note that the time scale for the response dynamics of each gene were set to be scaled to 1.
The FF-network type was more tolerant than the Direct type with simpler motifs. Indeed, the feed-forwaed motif has been reported to be tolerant to pulse-type inputs [2]. Interestingly, the Cooperative type was more tolerant to the pulse input than the FF-network type. In the Cooperative type, many ML-genes needed to change their expression levels in a sequential order before an output gene responded. An input gene stopped responding after a pulse interval, and the following ML-genes stopped responding before they stimulated the output gene. Thus, output genes avoided erroneous responses against a pulse input. Adaptive responses to unforeseen challenges All living organisms must survive unforeseen environmental changes to avoid extinction. However, evolution is not an easy solution for survival, as it requires the expense of time. If there are no individuals adapted to a new environment prior to beneficial evolutionary changes, they will go extinct. It is essential for cells to respond to a new environment to some degree without evolutions rewiring of the GRN, and this is referred to as an 'adaptive response to an unforeseen challenge (AUC)' [17].
We studied the possibility of the AUCs using the following procedure: we prepared a GRN model by doubling the number of input and output genes but maintained the number of ML-genes. Using only half of the input-output pairs, networks were evolved to optimize their fitness, that is, to give postulated output responses to the inputs. Here, half of the input and output genes were not used for the fitness, and hence they were not related to the original evolution. After the original evolution, we studied the responses of unused output genes when one of the unused input genes was activated. The response for this non-evolutionary input without further changes in the GRN is shown in Fig 5. Large responses to unforeseen inputs were observed in the Cooperative type networks, whereas there were almost no responses with the FF-network type. The Direct type had modest responses; however, this was because its core structure was so simple that even the addition of a random network could generate a response. In addition, such responses were rarer with larger network sizes. Hence, almost no response was expected in a network with many genes. The time course in the re-evolutionary process is shown in S1 Fig.

Resemblance of response dynamics between evolutionary and unforeseen inputs
One reason why networks of the Cooperative type realize the AUC is that many ML-genes exhibit similar dynamics in response to the activation of any input gene. They also behave in a similar manner even when an unforeseen input gene that is not used for original evolution is activated. Here, we quantitatively studied the resemblance of the response dynamics of ML-genes with different input signals using dynamic time warping (DTW), as shown in Fig 6 [32 -34]. DTW analyzes the similarity between the data of the two time-series taking into account their elongation and parallel shifts in the time axis direction.
For the Cooperative type, the ML-gene response dynamics to any input, including unforeseen ones, showed remarkable resemblance. These reflected cooperative behavior among the different inputs as well as among the different ML-genes. Although each input gene directly activated or inhibited a limited number of ML-genes as the network was sparse, almost all ML-genes subsequently responded through the entangled network with LEGI; some showed input specific responses, but most genes behaved independently of the input types. For the Direct and FF-network types, however, there were almost no such similarities among the different input signals that were observed, and only a small number of ML-genes exclusively responded to each input.

Dynamic mode decomposition of response dynamics
To investigate the origin of the resemblances of the response dynamics more closely, we analyzed the dynamics with dynamic mode decomposition (DMD) [35]. DMD decomposes time-series data with a large degree of freedom into a small number of dominant modes [36][37][38][39][40]. It corresponds to a principle component analysis (PCA) for temporal sequence data.
First, from a data series {x 0 , x 1 , · · · , x t }, we defined two time-series matrices X = (x 0 , x 1 , · · · , x t−1 ) and Y = (x 1 , x 2 , · · · , x t ) and assumed a relationship of Y = AX with an operator A. We calculated the eigenvalues and eigenvectors of A by utilizing the singular value decomposition. When the dimension of each time sample x i was D 0 , A was a square matrix of size D 0 and had D 0 pairs of eigenvalues and eigenvectors by definition. If only a small number (r) of singular values of X had much larger values than others that were close to 0, the singular value decomposition could be truncated to include only those r modes. The r modes dominantly influenced the dynam- ics of the data among the D 0 eigenvectors of A. These dominant eigenvectors corresponded to the primary axes in the PCA. How the number of dominant modes (r) differed according to the network type was investigated first (see Fig 7 (a)). For the Direct type, the number of dominant modes was almost the same as the number of respond-ing ML-genes with each input. For the FF-network type, most ML-genes showed no responses to the unforeseen inputs, and hence r = 1 truncation was just an artifact. In contrast, the Cooperative type showed responses that were represented by a smaller number of dominant modes. Its behavior could be described with fewer modes, although a larger number of ML-genes responded . r is defined as the number of singular values of a time-series matrix X that are larger than 1. Four different fitted networks were used for each parameter set, and 10 time-series data sets were computed for each network corresponding to each input, including both evolutionary and unforeseen inputs. (b) Distribution of DMD mode similarity between modes for an evolutionary input and modes for an unforeseen input (abscissa). After computing 5 dominant DMD modes (eigenvectors) for each input, the mode similarity was obtained as the average value of the inner products between eigenvectors for respective inputs. For the time-series data corresponding to 5 evolutionary inputs and 5 unforeseen inputs, that is, 5×5 = 25 combinations, the DMD mode similarity was computed to obtain the frequency histogram.
to each input. This result also indicated that ML-genes behaved cooperatively. It is of note that dimension reduction in high-dimensional adaptive states has been observed both in experiments and simulations [41,42], whereas the results of this investigation provided a reduction in the response dynamics.
Finally, the similarity in the dominant modes depending on the different inputs was investigated. With truncated eigenvectors from time-series data for each input, an inner product between eigenvectors for an evolutionary input and for an unforeseen input was calculated and the distribution of the products was computed. This inner product took a larger (smaller) value when the eigenvectors were aligned (orthogonal), that is, when the dominant modes were similar (diverse). As shown in Fig 7  (b), the Cooperative type showed the highest similarity, and the FF-network type had the lowest similarity. These results indicate that the response dynamics to an unforeseen input were like those that evolved with an evolutionary input for the Cooperative type. In other words, an unforeseen input was transmitted by utilizing an existing channel optimized for evolutionary inputs. Thus, the Cooperative type could respond to unforeseen inputs by activating any of the ML-genes in transfer of inputs.

Discussion
In this paper, we have compared the three types of GRNs that exhibit appropriate input-output responses: Direct, FF-network, and Cooperative types. It was demonstrated that the last type with complex entangled networks showed cooperative responses involving many genes and were robust to noise, errors in the inputs, and network changes. Furthermore, it also showed adaptive responses to unforeseen challenges. Hence, this type of network has biological advantages over the others. Even though the other two simple unidirectional circuits, either direct or feed-forward networks, could show rapid and accurate responses, the entangled network with cooperative responses could better deal with the fluctuating and unforeseen environmental changes.
It is of note that complex networks with cooperative responses emerge when the expression responses of genes are sloppy, in other words, have a low Hill coefficient. This contrasts with the unidirectional simple networks that emerge when the gene responses are sensitive and accurate (with high Hill coefficient). In general, even after evolution, the Hill coefficients are not very high, and the expression of each gene is sloppy. The Cooperative type networks evolved to compensate for the sloppy and unreliable response of each gene. Accurate responses by an ensemble of unreliable sloppy elements was originally proposed by von Neumann for the design of computers. Here, as a byproduct, the present gene expression dynamics were robust to perturbations and showed adaptive responses to unforeseen challenges [43].
The cellular networks are obviously entangled and composed of many elements that interact with each other. Cooperative responses have also been verified experimentally [44][45][46][47][48][49]. Gasch et al. studied the diverse environmental stress responses of yeast by transcriptome analysis. They reported that certain sets of genes (approximately 900 genes) exhibited similar responses to almost all environmental changes, while some genes showed unique response patterns to specific conditions only [44].
One of the most remarkable properties of cells is their ability to adapt to unforeseen changes in environmen-tal conditions [17][18][19][20]. Braun et al. found that yeast adapted to new environments that they had not previously experienced. This adaptive response occurred due to changes in their gene expression levels neither by means of prepared signaling networks nor by rewiring the networks via evolution. Despite the importance of such adaptive responses, how they are achieved remains elusive. The present investigation indicates a promising solution, by showing that the cooperative responses of the entangled networks induce similar responses for many genes for the experienced and unforeseen environmental inputs.
As a constructive experiment, Isalan et al. explored the effect of adding genes, i.e., new links to the GRN of E.coli [50,51]. The new gene could provide unforeseen inputs (or perturbations) to the existing GRN, and could provide similar gene expression patterns and outputs as those that were accounted for by the existing GRNs. In addition, the gene expression patterns fell into a small number of attractors. This behavior is consistent with the observations for the Cooperative type, exhibiting similarity in their responses to unforeseen challenges with that to preexisting inputs.
The origin of the correlated responses for many genes is due to local excitation specific to inputs and the global inhibition of many genes that are not specific to inputs. By considering the latter global response, the adaptive response to unforeseen challenges results. In biological systems, such local-excitation and global-inhibition (LEGI) is often observed, where inhibition occurs globally in space by using the global diffusion of inhibitors, for example, in chemotaxis pathways of eukaryotic cells [52] or Dictyostelium cells [53,54]. In this investigation, global inhibition took place not in real space but in network space. We expect that this LEGI behavior can be detected experimentally by global analysis of gene expression patterns and cellular pathways.
Cellular networks with many degrees of freedom are often studied by dividing them into small network motifs [55][56][57]. It is then assumed that each motif has a characteristic function, whereas the function of the whole network is given as a summation of the functions of the involved motifs. However, paths in the network are entangled, and motifs that work independently are difficult to extract. Moreover, it is not clear whether motifs embedded into an entangled network have the same function as in isolation [4,5].
At the same time, we do not necessarily state that the functions of network motifs are completely lost in an entangled network. Indeed, as already explained, the Cooperative type is robust against the changes in the parameter values because it includes the core structures of both the Direct and FF-network types (Fig 3). The Cooperative type can include such motifs that work in certain situations. In general studies, however, by taking advantage of the entangled network structure, the cooperative response is not decomposed into a combination of motifs.
In summary, the cooperative responses by entangled gene regulatory networks that are robust to noise and mutations can cope with unforeseen challenges, and this is essential for cell survival, as observed in recent microbial experiments.

Models
The gene regulatory network model is composed of N genes as nodes in a network, which are divided into three types: N I input genes receiving external inputs, N O output genes determining the fitness of the cell, and N M middle-layer genes (ML-genes) that transmit the input signals to the output genes (N I + N O + N M = N ). Through suitable normalization, the expression level of a gene is represented by a variable x i = [0, 1] (i = 1, . . . , N ), with the maximal expression level scaled to unity. The time evolution of each expression level is given as follows: where y i = I k δ ik + N j=1 C ij x j is the total signal received by the i th gene with δ ik as the Kronecker delta for k = 1, . . . , N I . I k shows the external input to the k th input gene. C ij represents the regulation from gene j to i, with 1 (excitatory), −1 (inhibitory), and 0 (nonexistent). y T denotes a constant threshold, and β determines the response sensitivity, which corresponds to the Hill coefficient. We assume that all genes in a network have the same y T and β values for simplicity.
For the evolution process, paths in the regulation matrix C ij were mutated, and C ij was selected according to the following fitness condition: for each given input, one predominant target among the output genes was expressed. Specifically, only the gene N −N O +k among the output genes N − N O + 1, · · · , N should respond upon the application of input I k (k = 1, . . . , N I ), where the response was given as the difference between the final (x f in i ) and initial (x ini i ) expression levels averaged over a time span, respectively, as With this fitness function, we evolved the network structure, that is the regulation matrix C ij by using a simple genetic algorithm. In the mutation process, we fixed the number of paths and swapped the connection C ij with a small mutation rate.
Supporting information S1 Fig. Time course in re-evolutionary process. Changes in fitness (ordinate) over a re-evolutionary course with an unforeseen input against the generation (abscissa). The networks are evolved to select a larger response to the new target gene corresponding to the unforeseen input. The averaged fitness from 125 different trials are plotted for each parameter set. Networks are re-evolved keeping the number of total paths from the fitted structure after the original evolution (a) or from its core structure (b).