Exploratory adaptation in large random networks

The contribution of outgoing hubs to the success of adaptation may reflect their ability to coordinate changes in a large set of affected nodes. In a network with a narrow distribution of out-degrees (without hubs), each node has the same relatively small influence as any other node. In the absence of a hierarchy in the extent of influence, irregular dynamic variation in the microscopic variables is unlikely to accumulate into a macroscopic coherent change in phenotype. On the other hand, the existence of a few hubs with a much broader influence can promote correlations between many downstream nodes, leading to an ability to drive a coherent change in a given direction. These effects may be related to other aspects of stability in network dynamics that vary with topology37,38,39.

Beyond the structural aspects promoting exploratory adaptation, the process of convergence appears to be complex and is characterized by an extremely broad distribution of times. Successful convergence likely depends on a delicate interplay between the space of possible network configurations, their connectivity properties and the typical timescales of their intrinsic dynamics.

While our model draws from neural network models40,41,42, it is substantially different in relying on purely stochastic exploration. In the language of learning theory, the ‘task' is modest: convergence to a stable attractor which satisfies a low-dimensional approximate constraint. Without exploration, this task could be fulfilled by chance with a very small probability. This probability increases dramatically by exploratory dynamics within a class of networks of a given structure. The ability to achieve high success rates without a need for complex computation or fine-tuning makes this type of adaptation particularly plausible for biological implementation. The relevance of similar processes in neural networks remains to be investigated.

Random network models were previously used to address evolutionary dynamics of gene regulation over many generations. These studies considered a population of networks undergoing random mutations and selection according to an assigned fitness21,43. In contrast, the model presented in the current study considers random variations over time within a single network, as an abstraction of a particular aspect of single cell adaptation within its lifetime. While these two approaches differ in timescales, level of organization and biological phenomena, it seems that they cannot be completely decoupled and that biological networks have basic properties that reflect on both contexts44. For example, in the context of selection in a population of networks, marked differences in evolutionary dynamics were found between homogeneous and SF networks45. In fact, the reproducible and exploratory responses in single cells, and the evolutionary processes at the population level, correspond to complementary aspects of gene-environment interactions at different scales3,46. A major future goal would be to integrate these aspects into a general picture of adaptive responses to diverse types of challenges over a broad range of timescales.

Interactions between the intracellular dynamical variables are governed by the network matrix W, defined as the element-wise (Hadamrd) product of the binary backbone, the adjacency matrix T, and a Gaussian random matrix J of connection strengths (equation (2)). We construct an ensemble of a given topology by sampling the connectivities of the backbone from a particular choice of in-degree and out-degree distributions, Pin(Kin) and Pout(Kout), and by sampling the random strengths of J independently from a Gaussian distribution. In practice, T is constructed first by randomly sampling a list of N out-going degrees from the distribution Pout(Kout) with ; and then sampling a list of N in-coming degrees from the distribution Pin(Kin) (again ), conditioned on the graphicality of the in- and out- degree sequences47. The network is then constructed from these sequences using the algorithm described in48.

Scale-free (SF) sequences are obtained by a discretization to the nearest integer of the continuous Pareto distribution P(K)=. Sampling SF degree sequences using the discrete Zeta distribution gives qualitatively similar results. Binomial sequences are drawn from a Binomial distribution , with p=. Exponential sequences are obtained by a discretization to the nearest integer of the continuous exponential distribution with . A Binomial degree sequence is implemented using MATLAB Binomial random number generator. Exponential and Scale-free sequences are implemented by a discretization of the continuous MATLAB Exponential and Generalized Pareto random number generators with parameters k=1/(γ−1), σ=a/(γ−1) and θ=a.

To compare adaptation performance between different ensembles, interaction matrices need to be properly normalized. In the study of uniform random matrices, the elements are usually normalized such that their variance is , providing a well-defined thermodynamic limit N→∞ in which the matrix eigenvalues of are uniformly distributed within a disc of size g0 in the complex plane49,50.

In our model, the initial interaction matrix is defined as a random Gaussian matrix with mean 0 and variance , being the average connectivity. Neglecting correlations in the adjacency matrix T, the variance of its elements is , which implies Var(Wij)≈. In principle both finite-size effects and correlations in Wij result in deviations from a uniform distribution of eigenvalues in the circle. However empirically we find that for matrices of relevant size, the spectral radius of W is still ∼g0, establishing a basis for comparison between the different ensembles based on spectral radius. We note however that the eigenvalue distribution is far from being uniform (see Supplementary Note 1, Empirical spectrum of interaction matrices W).

Another model component that needs to be normalized for proper comparison is the macroscopic phenotype y(x)=b·x. The arbitrary weight vector b is characterized by a degree of sparseness c, the fraction of non-zero components, ; and by the typical magnitude of those components. In order to compare between networks of different size and weight vectors of different sparseness, the variance of the non-zero components is scaled by their number, cN and by the matrix gain g02. The non-zero components of b are thus distributed as , where α is a single parameter determining the scale of phenotype fluctuations in different network sizes and gains (See Supplementary Note 1, Distributions of phenotype y).

Convergence fractions were computed over 2,000 time steps in samples of 500 networks drawn from specified in- and out-degree distributions, averaging over T, J0 and x0. For fully or sparsely connected homogeneous random networks of size N=1,500, the CF is close to zero (not shown). Alternative ensemble definitions (for example, keeping T fixed) do not change the main results (see Supplementary Note 1, Convergence of different network ensembles).

The saturating function is defined as an element-wise function φ(xj)=tanh(xj) operating separately on each of the components of x. Model results are insensitive to the exact shape of this function (Supplementary Note 2, Robustness of model to saturating function φ) and to placing the saturation inside or outside of the interactions (Supplementary Note 2, Robustness of model to position of saturating function φ).

The mismatch function is defined here as , a symmetric sigmoid around y*, where ɛ=3 controls the size of the low-mismatch ‘comfort-zone' around y*, μ=0.01 the steepness of the sigmoid and 0=2 its maximal value. Main model results are insensitive to the exact shape of this function as long as it has a flat region with zero or very low mismatch around y*. (see Supplementary Note 2, Robustness of model to mismatch function ).

The data that support the findings of this study are available from the corresponding author upon reasonable request.

How to cite this article: Schreier, H. I. et al. Exploratory adaptation in large random networks. Nat. Commun. 8, 14826 doi: 10.1038/ncomms14826 (2017).

Publisher's note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

We thank O. Barak, E. Braun, R. Meir and M. Stern for valuable discussions and S. Marom, A. Rivkind, L. Geyrhofer and H. Keren for critical reading of the manuscript.

The data that support the findings of this study are available from the corresponding author upon reasonable request.