Exploratory adaptation in large random networks
The capacity of cells and organisms to respond to challenging conditions in a repeatable manner is limited by a finite repertoire of pre-evolved adaptive responses. Beyond this capacity, cells can use exploratory dynamics to cope with a much broader array of conditions. However, the process of adaptation by exploratory dynamics within the lifetime of a cell is not well understood. Here we demonstrate the feasibility of exploratory adaptation in a high-dimensional network model of gene regulation. Exploration is initiated by failure to comply with a constraint and is implemented by random sampling of network configurations. It ceases if and when the network reaches a stable state satisfying the constraint. We find that successful convergence (adaptation) in high dimensions requires outgoing network hubs and is enhanced by their auto-regulation. The ability of these empirically validated features of gene regulatory networks to support exploratory adaptation without fine-tuning, makes it plausible for biological implementation.
The ability to organize a large number of interacting processes into persistently viable states in a dynamic environment is a striking property of cells and organisms. Many frequently encountered perturbations (temperature, osmotic pressure, starvation and more), trigger reproducible adaptive responses1,2,3. These were assimilated into the organism by variation and selection over evolutionary time. Despite the large number and flexible nature of these responses, they span a finite repertoire of actions and cannot address all possible scenarios of novel conditions. Indeed, cells may encounter severe, unforeseen situations within their lifetime, for which no effective response is available. To survive such challenges, a different type of ad-hoc response can be employed, utilizing exploratory dynamics4,5,6,7,8.
The capacity to withstand unforeseen conditions was recently demonstrated and studied using dedicated experimental models of novel challenge in yeast9,10,11,12 and flies13. Adaptive responses exposed in these experiments involved transient changes in the expression of hundreds of genes, followed by convergence to altered patterns of expression. Analysis of repeated experiments showed that a large fraction of the transcriptional response can vary substantially across replicate trajectories of adaptation10,12. These findings suggest that coping with unforeseen challenges within one or a few generations relies on induction of exploratory changes in gene regulation over time in an individual5,6.
Several properties of gene regulatory networks may support such exploratory adaptation. These include a large number of potential interactions between genes14, context-dependent plasticity of interactions15,16,17,18 and multiplicity of microscopic configurations consistent with a given phenotype19. Despite these properties, the feasibility of acquiring adaptive phenotypes by random exploration within a single organism remains speculative and poorly understood. In particular, it is not known how exploration may converge rapidly enough in the high dimensional space of possible configurations? what determines the efficiency of this exploration? and what ensures the stabilization of new phenotpes?
Here we address these open questions by introducing a network model of gene regulation, which demonstrates the capacity to adapt by exploratory dynamics in a single cell (as opposed to selection on existing variation in a population). Exploration is triggered by failure to satisfy a newly-imposed external demand, and is implemented by a random walk in the space of network configurations. Exploration relaxes if and when the system reaches a stable state satisfying this demand. We show that the success of this exploratory adaptation in high dimension requires that the network include outgoing hubs. Adaptive capability is further enhanced by autregulation of these outgoing hubs. Since these are both well-known properties of gene regulatory networks, our findings establish a basis for a biologically plausible mode of adaptation by exploratory dynamics.
To investigate the feasibility of exploratory adaptation, we introduce a model of gene regulatory dynamics incorporating random changes over time in a single network. The model consists of a large number, N, of microscopic components x=(x1, x2 ... xN), governed by the following nonlinear equation of motion (Fig. 1a):
where W is a random matrix, representing the intracellular network of interactions; φ(x) an element-wise saturating function restricting the dynamic range of the variables; and the relaxation rates are set to unity. Previous work has used similar equations to address evolutionary aspects of gene regulation20,21 as well as interactions and relaxation in neuronal networks22. Most studies have focused on networks with uniform (full or sparse) connectivity; much less is known about the dynamics for networks with non-uniform topological structures, which may be of relevance to gene regulation.
Here we consider sparse random networks with different types of topological properties. For all cases, the interaction matrix W is composed of an element-wise (Hadamard) product,
where T is a random topological backbone (adjacency) matrix with binary (0/1) entries representing potential interactions between network elements; and J is a random matrix specifying the actual interaction strengths. To represent context-dependent regulatory plasticity, we assume that the backbone remains fixed, whereas the interaction strengths are plastic and amenable to change over time. We will emphasize below network sizes and topological structures that are relevant to gene regulatory networks.
On a macroscopic level we consider a cellular phenotype, y, which depends on the microscopic components and can affect the cell's functionality and state of stress. We define this phenotype as a linear combination of microscopic variables
with an arbitrary vector of coefficients b. To model an unforeseen challenge, the system is subjected to an arbitrary contstraint of maintaining the phenotype in a given range y(t)≈y*. Importantly, any given value of the phenotype can be realized by many alternative microscopic combinations.
Deviation from compliance with the constraint is represented by a global cellular function (y−y), corresponding to the level of mismatch between the current phenotype and the demand. This mismatch is effectively zero inside a ‘comfort zone' of size ɛ around y and increases sharply beyond it. Biologically, the comfort zone can be interperted as a range of phenotypes that can be tolerated in a given environment without invoking significant stress. This is represented mathematically by a range of values which satisfy the constraint (in contrast to many optimization problems which require adherence to a specific value).
When the phenotype deviates from the comfort zone, the mismatch drives an exploratory search, realized by small random changes in the interaction strengths, forming a random walk in the elements of the matrix J:
where is the standard Wiener process. The amplitude of the random walk is controlled by a scale parameter, D, and the mismatch level, . These random changes can arise from diverse sources of variation affecting the levels of transcription regulators3,23,24, as well as regulatory interactions (for example, alternative splicing, conformations of transcription factors and their post-translational modifications17,18).
The random walk constitutes an exploratory search for network configurations in which the dynamical system in equation (1) satisfies the constraint in a stable manner. Random occurrence of such a configuration decreases the search amplitude, thereby promoting relaxation by reducing the drive for exploration6,25. Convergence of this process to a stable state satisfying the constraint is not a priori guaranteed. Intuitively, it may be expected that randomly varying a large number of parameters in a nonlinear high-dimensional system will cause the dynamics to diverge. Surprisingly, we find that the adaptation process can in fact converge; however, as shown below, convergence depends on key properties of the network.
An example of adaptive convergence is shown in Fig. 1b–d. At t=0, the system is confronted with a demand and starts an exploratory process in which the connection strengths are slowly modified. Figure 1b displays the time trajectories of four of these connection strengths. During this exploration, the microscopic variables, x, and the phenotype, y, exhibit highly irregular behavior, rapidly sampling a large dynamic range (Fig. 1c,d respectively). At t∼400, the system manages to stably reduce the mismatch to zero and converges to a fixed point (Fig. 1). In some cases the dynamics converges to a small-amplitude limit-cycle (Supplementary Note 3, Convergence to a limit cycle), and remain within the comfort zone ±ɛ around y. The state of convergence is found to be a stable attractor that is robust against small perturbations of the dynamic variables, x, and the interactions strengths, Wij (Supplementary Note 3, Stability of the adapted state). The differences between the amplitude of temporal changes in Fig. 1b–d reflects the separation of timescales between the slowly accumulating changes in interaction strengths, governed by the small value of D in equation (4), and the intrinsic dynamics of equation (1).
To investigate the dependence of exploratory adaptation on network topology we constructed random matrix ensembles with different topological backbones, manifested by distinct in- and out-going degree distributions26 (detailed in the Methods section). Each ensemble was evaluated with respect to the probability of convergence, estimated as the fraction of simulations which converged within a given time window. Figure 2a compares ensembles of networks with in- and out-degrees drawn from Binomial (Binom), Exponential (Exp) and Scale-Free (SF) distributions. It shows high fractions of convergence, 0.5 or higher, only for ensembles with SF out-degree distributions. In contrast, the in-degree distribution affects convergence only mildly. For example, the convergence fraction (CF) of networks with SF out-degree and Binomial in-degree distributions (dark blue) is 0.5, whereas it is only 0.03 in the transposed case (light blue). This asymmetry between outgoing and incoming connections indicates that convergence of exploratory adaptation does not rely on spectral properties of the interaction matrix ensemble.
Analysis of convergence as a function of network size shows that the effect of topology becomes pronounced for large networks (Fig. 2b). The CF in small to intermediate-sized networks (N≲200) is higher and relatively independent of topology. However, as N increases towards sizes that are relevant to genetic networks, the benefit of having SF out-degree distribution becomes progressively prominent.
Among the topological ensembles tested, an outgoing SF degree distribution was found to be crucial for convergence in large enough networks. Such distributions are characterized by a broad range of heterogeneous connectivities, with a small number of extremely highly connected nodes (hubs). To evaluate the relative contribution of outgoing hubs to convergence within this ensemble, we ranked the backbones of the connectivity matrices drawn from the SF-Binom distributions according to the out-degree of the largest hub. Figure 3a shows that the CF increases with the connectivity of the largest outgoing hub. As a second approach to characterize hub contribution, we deleted a small number of outgoing hubs from these networks27; this leads to a significant reduction in CF that is not observed upon removal of randomly chosen nodes (Fig. 3b).
These results indicate that, in networks from the SF-Binom ensemble, outgoing hubs have a major positive influence on the success of exploration. We therefore asked whether the addition of a few hubs to an otherwise poorly converging ensemble is enough to induce significant convergence. Figure 3c indeed shows that addition of as few as 8 hubs to a Binom-Binom ensemble increases the CF from zero to about 0.4.
These findings are in-line with reported properties of gene regulatory networks, particularly the existence of ‘master regulatory' transcription factors that control the expression of hundreds of other genes28,29,30. Since many of these master regulators are also autoregulated31, we evaluated the influence of hub autoregulation on the success of exploratory adaptation in our model. Figure 3d shows that autoregulation of the leading hubs in the SF-Binom ensemble further increases the CFs.
Since autoregulation motifs are commonly observed in gene regulatory networks (not only in hubs)32, we investigated whether these motifs could also contribute to convergence when over-represented uniformly throughout the network. Figure 4 depicts the results of adding such motifs randomly to 10% of the nodes in networks from the SF-Binom and Binom-SF ensembles. It is seen that positive autoregulation enhances convergence for intermediate sized networks (N=1,000) in both ensembles; this effect is particularly notable for the Binom-SF ensemble, which has small CF without these motifs. This contribution, however, decreases with network size and vanishes in the same type of networks with N=3,000. We conclude that the presence of autoregulatory motifs ranodmly positioned in the network cannot substitute for hub contribution in the limit of very large networks. These results highlight the interplay of several networks properties in exploratory adaptation: network size, topology and autoregulatory motifs. The addition of common network motifs other than autoregulation did not lead to a conclusive effect on convergence (Supplementary Note 3, Dependence of convergence on network motifs).
We investigated how the capacity to adapt is affected by various model parameters. To examine the dependence on the severity of the constraint, we varied the size of the comfort zone ɛ. Figure 5a reveals a sharp decrease of the CF as ɛ is reduced, indicating that a non-vanishing comfort zone is crucial for successful exploratory adaptation. This requirement is biologically plausible, as one expects a range of phenotypes capable of accommodating a given environment rather than a unique optimal phenotype. Another way of increasing the adaptation challenge is by shifting the required phenotypic range away from the origin. Reaching a shifted region is challenging because it is more rarely visited by spontaneous dynamics (Fig. 5b, grey curve). Figure 5b indeed shows that the CF decreases as y* moves away from zero (blue curve). Importantly however, it remains much larger than the probability of encountering the required phenotype spontaneously. For example, a non-negligible convergence (CF∼0.2) is observed even for an interval around which is spontaneously encountered with probability of 0.02.
To evaluate the sensitivity of adaptation to exploration speed, we varied the effective diffusion coefficient in the space of connection strengths, D. Figure 5c shows that a non-zero convergence fraction is achieved for a wide range of this parameter and remains between 0.2 and 0.7 over more than 5 orders of magnitude. As the value of D increases beyond a certain level where the separation of timescales ceases to hold, the convergence fraction decreases rapidly.
For a given adjacency matrix T, interactions within the network are determined by the connections strengths, Jij. These are initially drawn from a Gaussian distribution with a zero mean and a given s.d. The s.d. normalized to network size, g0 (also called network gain; for details see Methods) determines the contribution of the first versus second term in equation (1). In large homogeneous networks, this parameter has a strong effect on the dynamics of equation (1) (ref. 33). In contrast, we find that the capacity to adapt by exploration in our model is relatively weakly dependent on g0 (Fig. 5d).
The analysis presented so far was based on convergence fractions within a fixed time interval. To characterize the temporal aspects of exploratory adaptation, we evaluated the distribution of convergence times in repeated simulations. Figure 6 reveals a broad distribution (CV≈1.1), well fitted by a stretched exponential (see Supplementary Note 3, Stretched exponential fit to the distribution of convergence times). Such distributions are common in complex systems34 and were suggested to reflect a hierarchy of timescales35. While the general shape of the distributions were similar in all topological ensembles tested, networks with SF out-degree distributions typically converged faster than their transposed counterparts (Fig. 6a). Moreover, deletion of a small number of leading outgoing hubs causes a significant shift towards longer convergence times (Fig. 6b). Thus, networks with larger heterogeneity in out-degrees are both more likely to converge within a given time window (Figs 2 and 3), and typically converge faster (Fig. 6).
In the typical example shown in Fig. 1, exploratory dynamics culminates in reduction of drive for exploration and convergence to a stable attractor of equation (1). The significant differences between adaptive performance of network ensembles (Fig. 2a,b) may reflect the abundance of networks supporting relaxation to attractors in the different ensembles. Previous work has shown that for networks with uniform degree distributions and sufficiently strong interactions, the number of attractors of equation (1) decreases with network size and vanishes in the limit of infinite size (leading to chaotic motion only33). A related result was recently found for Boolean networks36. It is not known, however, how the number of attractors scales with system size for networks of arbitrary topological structure.
To address this question, we simulated many independent networks in each ensemble and estimated the fraction which relaxed to fixed points without exploration or feedback (equation (1) alone). For any given network, the probability of relaxation to a fixed point was largely insensitive to the initial conditions in x-space (not shown). With that in mind we computed, for each topological ensemble, the fraction of networks supporting relaxation within a given time window, starting with random initial conditions. This measure is analogous to the CF used in Fig. 2, but without any constraint, feedback or random walk in connection strengths. To highlight the dependence on network size we extended the simulations up to N=10,000. Figure 7a reveals topology-dependent differences that are qualitatively in line with the ability for exploratory adaption shown above (Fig. 2b). This suggests that a substantial contribution to successful adaptation is indeed provided by a high abundance of networks exhibiting fixed points in their dynamics.
For each network ensemble that supports fixed points, we further analysed the distribution of relaxation times into these fixed points. Figure 7b demonstrates the effect of topology by comparing the SF-Exp ensemble to the transposed Exp-SF. It shows that networks with SF-out degree distribution typically support faster relaxation to their respective fixed points. This may allow the adaptation to converge before exploration has had a chance to significantly modify network connections. Further work is required to test this hypothesis and to broaden the theoretical understanding of these dynamics in random ensembles with heterogeneous topologies.
Overall, we have introduced a model of exploratory adaptation driven by mismatch between an internal global variable and an external constraint. Adaptation is achieved by a purely exploratory process which relies on the plasticity of regulatory interactions17,18. Our model was described in terms of gene regulation but could equally well represent adaptation in other cellular interactions, such as the protein-protein interaction network. We have found that convergence of exploratory adaptation depends crucially on structural properties of the network. It requires the existence of outgoing hubs and is enhanced by auto-regulation of these hubs. These results offer an important, but hitherto unrealized, rationale for the overwhelming abundance of autoregulation motifs on master regulatory transcription factors31. These master regulators act as network hubs by virtue of the large numbers of their downstream gene targets. Our findings show that autoregulation of such hubs improves their ability to drive the network into a stable state which satisfies a phenotypic demand.