Explorability and the origin of network sparsity in living systems

The increasing volume of ecologically and biologically relevant data has revealed a wide collection of emergent patterns in living systems. Analysing different data sets, ranging from metabolic gene-regulatory to species interaction networks, we find that these networks are sparse, i.e. the percentage of the active interactions scales inversely proportional to the system size. To explain the origin of this puzzling common characteristic, we introduce the new concept of explorability: a measure of the ability of an interacting system to adapt to newly intervening changes. We show that sparsity is an emergent property resulting from optimising both explorability and dynamical robustness, i.e. the capacity of the system to remain stable after perturbations of the underlying dynamics. Networks with higher connectivities lead to an incremental difficulty to find better values for both the explorability and dynamical robustness, associated with the fine-tuning of the newly added interactions. A relevant characteristic of our solution is its scale invariance, i.e., it remains optimal when several communities are assembled together. Connectivity is also a key ingredient in determining ecosystem stability and our proposed solution contributes to solving May’s celebrated complexity-stability paradox.

In inanimate matter, elementary units, such as spins or particles, always have their mutual interactions turned on (with the intensity decaying with their relative distance). Thus, the interaction network is dense, with all connections present, i.e. particles do not have the freedom to adjust or change their interactions unless they change their relative distances. In contrast, living systems are composed of interacting entities, such as genes [1][2][3][4] , metabolites 1,5,6 , individuals [7][8][9] , and species 4,10-14 , with the ability to rearrange and tune their own interactions in order to achieve a desired output 1 . Indeed, thanks to advances in experimental techniques, which are generating an increasing volume of publicly available ecologically and biologically relevant data, several studies indicate that interaction networks in living systems possess a non-random architecture characterised by the emergence of recurrent patterns and regularities 10,11,15,16 .
Analysing different data sets of ecological, gene-regulatory, metabolic and other biological interaction networks 1,2,11,12,[17][18][19] , (see Supplementary Information for details), we find that one ubiquitous emergent pattern is sparsity, i.e. the percentage of the active interactions (connectivity) scales inversely proportional to the system size (illustrated in Fig. 1). For example, in the case of ecological systems, species interact selectively even when they coexist at short distances and most of the interactions are turned off. A generic system formed by S interacting units may have a maximum number of interactions equal to S 2 (including self-interactions), i.e. a connectivity C (defined as the fraction of active interactions) equal to 1. On the other hand, the minimum number of interactions that guarantees that the interaction network is connected is of the order of S, that is ∼ C S 1/ , corresponding to the percolation threshold of random networks 20 . Thus, in this range of possible connectivities, it is quite surprising that the observed ones in the analysed interaction networks all correspond to the lowest possible values. However, it is not known if this recurrent property gives any advantage or reward to the system and a theoretical framework to understand the origin of sparsity is still lacking.
Guided by such an intriguing observation, the main goal of our work is to shed some light on why this pattern emerges and to study, from a theoretical point of view, if the sparsity of interactions confers any advantage to the system. In this context, variational principles have been proven to be a useful tool to elucidate some of the recurrent patterns in nature 11,21 , In the same vein, in this work we propose an optimisation approach to describe the role of active interactions in living systems. We show that sparse networks offer, at the same time, a maximum capability of the system to visit as many stable attractors as possible by simply tuning the interaction strengths (explorability), as well as the largest robustness of the underlying dynamics, guaranteeing that such attractors remain stable (dynamical robustness).

Results
Mathematical framework. We consider a system composed of S nodes (e.g. species, metabolites, genes) characterised by dynamical variables, = .
.. (e.g. populations, concentrations, levels of expression), following a generalised Lotka-Volterra (GLV) dynamics: In the simplest case , we recover the classic S species Lotka-Volterra equations and we refer to the parameter α i as the growth rate. In this case the interaction of node i with node j is encoded in the matrix element w ij , whose diagonal entries set the scale of the interaction strengths, which for the sake of simplicity 22 we set to −1. For convenience, we also introduce the adjacency matrix A ij , whose entries are 1 if the corresponding ≠ w 0 ij and 0 otherwise. A non-trivial stationary point of the dynamics of Eq. (1), ⁎ x , is determined by the interactions within the system, i.e. when ∑ = w 1 , and its stability is guaranteed if all the eigenvalues of the Jacobian matrix evaluated at this point, = ⁎ J x w ij i ij , have a negative real part. GLV dynamics have been used to model the time evolution of ecological systems 22,23 , human microbiota dynamics 24 , gene expression 25 and other biological systems 26 , where x i represents the density of the i-th species, and therefore we focus on the stable and feasible stationary solutions of the dynamics 19,22 For simplicity, here we focus on the simpler case with , but in the Supplementary Information we recast all the results for the more general non-linear dynamics, Eq. (1).
Explorability. Depending on the model parameters, the dynamics of (1) exhibits different fixed points, , that may correspond to feasible/non-feasible and either stable/unstable solutions. We do not study more complicated possibilities such as limit cycles or chaotic strange attractors, as this would require to go beyond linear stability analysis. Focusing on the interaction strengths as our tuning parameters, by varying them we can pass from one fixed point to another, and change its feasibility and stability. In brief, we take the explorability as the volume of feasible and stable fixed points spanned by modifying the link weights, while all other model parameters are kept fixed. The explorability resembles the concepts of robustness and adaptability in the context of evolutionary dynamics, which make up the number of 'phenotypes' (attractors of the dynamics) that can be reached by mutations in the space of 'genotypes' (interaction networks) 27,28 , Nevertheless, in our case the system has only one stable and feasible fixed point ('phenotype') once the weights w ij are fixed (w is invertible). Therefore, for any change of a link weight (a 'mutation' of w) we have a different fixed point ('phenotype'), i.e. in our setting there are not 'neutral' mutations.
More specifically, we can define the explorability for interacting systems described by a GLV dynamics Eq. (1) and for a given topology (i.e., an ensemble of interaction networks having the same adjacency matrix A, but in general different links weights w ij ) of the interaction network. The explorability V E is the volume in  S spanned by all the feasible and stable fixed points as one varies w ij while keeping all other parameters fixed (see Fig. 2). Notice that, with such a definition, the fully connected network has the largest possible explorability, since any other topology is attainable by making some of the matrix entries arbitrarily close to zero. However, might the optimal or quasi-optimal solutions indeed be the ones where most of the interactions are turned off, as suggested by the observational data? Moreover, in the fully connected case, many interaction parameters have to be specified (there are S 2 matrix elements that can be varied), and then spanning all the possible fixed points becomes a complicated, fine-tuning problem, which does not seem to be feasible in biological systems 29 . Therefore, we pose the following questions: what is the relationship between explorability and the interaction network topology? Is there an optimal network structure that maximises explorability? To answer these questions, we started by analysing the extreme case of a sparse topology with just S links, i.e. a tree with one loop with connectivity C = 2/S (see Fig. 2) (the factor 2 comes from the fact that we also count the self-interactions). For the sake of simplicity in what follows we will refer to this topology as the tree-like network. Even in this simple case, measuring the explorability requires us to scrutinise an S-dimensional space of parameters (corresponding to the S entries of the interaction matrix). Furthermore, one still has to choose the values of α i , which are intrinsic parameters of the dynamics (e.g. species growth rates) and should be set a priori (in contrast to the interactions w ij , which are tuning parameters in our approach). For this reason, we introduced various degrees of approximations in our setting. We first considered the simplest uni-parametric case, α α = i , that we can fix without a loss of generality to α = 1 (we checked that our conclusions are still valid for other choices of α i , see below). In addition, we first restricted the analysis to the subspace of fixed points with homogeneous components, i.e. = ⁎ ⁎ x x i . Under these approximations, computing the explorability becomes a much simpler task and we were able to develop an analytical solution to this problem (see Methods). This approach, although a priori seems too drastic, leads to rather reasonable estimates of the explorability. As a second step, we enlarged the region explored by introducing some heterogeneity into the components of x*.
Most probably, the volume of feasible and stable fixed points may become infinite. However, we are always interested in comparing such volumes for different topologies. Indeed, in the simple homogeneous situation, we observed Measuring explorability and dynamical robustness. We considered a system of S nodes representing species, genes, metabolites,… (S = 7 in the picture), whose state = . ..
, with w ij encoding the network structure and the strength of the interactions. (Panel A) We started from a tree-like network, i.e. a tree with one loop, with S links, represented by the bluecolored links, and, for such a topology, we searched for the feasible and stable fixed points as we varied the interaction strengths. (Panel B) The spanned volume sketched by the blue-shaded region (a 2D projection of the S-dimensional space) corresponds to the network explorability. (Panel C) The dynamical stability quantifies the robustness of the system to perturbations of the dynamics itself, ) ( ) , and should not be confused with the standard stability or resilience that measures the response to perturbations on the fixed point (which is already implicit in the measure of explorability). Starting from an interaction matrix whose fixed point of the underlying dynamics is at the edge of stability, denoted by w ij edge (for which λ = R( ) 0 max , black dot in the picture), the dynamics was perturbed, and the stability of the new fixed point was evaluated in order to test the dynamical robustness of the system. As demonstrated in the Supplementary Information, this can be simply encompassed by computing the principal eigenvalue λ′ R( ) max of the perturbed Jacobian matrix, ξ ′ = ′ J w ij i ij , where ξ′ depends on the details of the perturbation, that we take to be a random vector. The histogram of λ′ R( ) max is sketched in panel C as the dynamics ξ′ is varied. Finally, we increased the connectivity of the network by including additional fixed strengths, ε ij (red edges in the graph), to the network. The same previous analysis was performed again by varying the strengths of the blue links and the corresponding results are shown in red in the panels on the right. By randomly sampling the location and strengths of the added links, we investigate if the explorability and dynamical robustness statistically increase or diminish for higher values of the connectivity. that in almost all the cases (100 per cent for tree-like networks and, e.g., more than 98 per cent for C = 0.5), we can identify two regions along = ⁎ ⁎ x x i : fixed points become unstable for small x* and stable for large x*, with a marginally stable fixed point intersecting at a single value ⁎ x c , for which λ λ = = = ..
, where λ i are the eigenvalues of the Jacobian matrix J (see Supplementary Information for more details). Therefore, we can take as a proxy for the explorability, where V 0 is a sufficiently large constant (V 0 = 1 in our analysis), which allows for comparison between different topologies. With this definition, we can compute analytically the explorability of a tree-like network with S links (see Supplementary Information), finding that, among all the possible tree-like topologies, the one with just a loop composed of three nodes leads to the optimal explorability, = V 2/3 E . This structure (which we refer to as the optimal tree-like network) constitutes our reference network when increasing the connectivity.
As a third step, we analysed the explorability of networks with higher connectivities. This enormously increases the number of matrix entries that have to be modified when computing the spanned volume of feasible and stable fixed points. For this reason, we adopted the following approach (see Fig. 2): starting from the tree-like topology, we introduced additional links to the tree-like topology of weights ε ij for any extra links between nodes i and j, and then computed the explorability, fixing the values of all ε ij and tuning the other matrix elements.
Sampling different values of the added link weights (but not their locations), we can construct a histogram of the explorabilities, In addition, we are not interested in distinguishing different topologies with the same connectivity C, so we also sampled over different locations of the added links, leading to | P V C ( ) E (see Methods for technical details). For numerical reasons, we have tested that our results are robust for network sizes S < 100. However, it is important to mention that network size is a significant factor in the spectrum of fixed points 30 and that more complex phenomenology could be found in very large systems.
Numerical results are represented in Fig. 3 (top left panel), illustrating that the explorability of the optimal tree-like network is indeed statistically higher than the one for denser networks. Furthermore, the average explorability decreases as the connectivity of the interaction network increases (Fig. 3, top right panel). We were able to prove this result for some particular topologies of small networks, for which the explorability can be calculated analytically (see Supplementary Information). In conclusion, our results suggest that, on average, explorability decays with the connectivity of the system, and therefore, sparse networks generally lead to higher values of explorability.
Dynamical robustness. Another crucial property of complex interacting systems is their robustness to perturbations 31,32 , Understanding the role of network architecture in the stability of a system with many degrees of freedom is an important challenge, since it impacts on our capacity both to prevent system failures and to design more robust networks to tolerate perturbations to the system dynamics.
The standard measure of stability (known as asymptotic resilience in ecology 33,34 ,) is defined as the capacity of the system to return to the original stationary state after a perturbation of it, δ + ⁎ x x, while the dynamics is kept fixed. Let us note that our definition of explorability already takes into account this kind of stability, i.e. explorability is defined only for stable system dynamics.
Alternatively, we can study how the stability of the system is modified as a result of a perturbed dynamics, where δG and δF represent the perturbations with respect to the original dynamics. This can be understood as including further non-linear effects that were not present before. As a consequence of this kind of perturbation, both the original stationary states and their degree of stability are modified. We then quantify the capacity of the system to re-organize after a perturbation of the dynamics such that the new stationary state of the system is close to the original one and still stable. We refer to this as the dynamical robustness of the system (to avoid confusing this new measure of stability with the standard resilience).
We found the pleasing result that the Jacobian matrix evaluated at the new stationary state, ′ J ij , retained a similar form than for the original dynamics, i.e. ξ ′ = ′ J w ij i ij , where ξ′ depends on the specific details of the perturbed dynamics (see Supplementary Information). Focusing on the worst case, which corresponds to marginally stable fixed points for which λ = R( ) 0 max (denoted by w edge ), we perturbed the dynamics and computed the maximum real part of the eigenvalues of ′ J ij at the new stationary point, λ′ R( ) max . For a given deterministic perturbation of the dynamics ξ′ and topology, λ′ R( ) max can be taken as a measure of the dynamical robustness (see Fig. 2, panel C). Since we wanted to keep the analysis as general as possible, we studied the dynamical robustness against random perturbations of the dynamics (generated from a distribution ξ′ P( )). In particular, for each connectivity C, we fixed the additional links to ε ij and then looked for the set of matrices ε w ({ }) ij edge at the edge of stability (for which, by definition, λ = R( ) 0 max ). Taking different realisations of ε ij and ξ′ i , we compared the distributions λ′ | R P C ( ( ) ) max (see Fig. 2). We can then define a statistical measure R of the dynamical robustness by taking the value of λ − ′ R( ) max located at the fifth percentile, 5 . In this way, we have an indicator of how much stability could be gained under random perturbations of the dynamics. Other choices of this measure can be taken, e.g. based on the 10th, 20th and 50th percentiles, leading qualitatively to the same conclusion. However, lower percentile values generally enhance the differences between topologies (see Supplementary Information).
The bottom left panel of Fig. 3 shows the histogram of λ′ R( ) max for different connectivities. Again, it can be seen that the case of a tree-like network leads to the best performance, whereby the fixed point of the perturbed dynamics generally becomes stable, i.e. λ′ R( ) max is negative, and the modulus reaches larger values than in the corresponding case of networks with higher connectivity. Therefore, our analysis shows that sparse tree-like networks have both a larger explorability and a larger dynamical robustness than random networks with higher connectivity.
Optimisation approach. We then went one step further and compared the explorability and dynamical robustness of tree-like networks with graphs constructed via an optimisation process rather than randomly generated.
In the optimisation, weight values ε ij were changed accordingly to a stochastic hill climbing algorithm 35 , whereas their locations were kept fixed. Introducing a random Gaussian perturbation with zero mean and standard deviation δε in all ε ij , the new configuration was accepted if the quantity to optimise (either V E or R) increased, and the process was iterated for T time steps. Our results are robust for different choices of δε in the range [10 −3 , 10 −1 ]. We restricted our analysis to the homogeneous case ( ), using network sizes of S = 10 to facilitate the convergence of the optimisation algorithms. Still, the landscape of ε V ( ) E ij and ε R( ) ij appeared to be highly irregular with many local minima when increasing the connectivity. The corresponding eigenvalue with the largest real part ( λ′ R( ) max ) gave the degree of stability for the new fixed point, from which we measured the dynamical robustness of the system (see Methods). We generated 10 3 configurations of ε ij 's as done for the top panel (ε = 0 ij corresponds to a tree-like network), and for each one, 10 3 values of ξ′ (encoding the perturbation of the dynamics) from a uniform distribution in [0,1], and then computed the system response through the distribution λ′ . Both homogeneous and nonhomogeneous settings were analyzed for the tree-like network and = .
C 0 5 (same color code as for the top left panel). In all cases, we found that increasing the network connectivity shifted the distribution of R λ ′ ( ) max towards the less stable region. Qualitatively similar results were obtained for larger values of σ ε , σ x and σ α (see Supplementary Information). (Bottom right) Dynamical robustness R (estimated as the value of λ′ R( ) max located at the fifth percentile of the distribution,  Figure 4 represents the initial random-generated networks with connectivities in the range ∈ .
. C [0 25, 0 75] (blue dots) and the corresponding optimised values for the explorability and the dynamical robustness (green and magenta dots, respectively). Qualitatively similar results are found using different schedules of the simulated annealing algorithm 36 . The explorability reached for networks with connectivities larger than the one for tree-like networks turned out to be very close to the corresponding value for the optimal tree-like topology. However, in general, such networks exhibited low values of dynamical robustness. Similarly, when optimising the dynamical robustness we ended up with values very close to that of the optimal tree-like network, but remarkably, without improving the explorability. Indeed, in all the cases that we have analyzed, during the optimization of one property, the other one either remained unaltered or decreased. However, for denser networks the optimization becomes harder due to the presence of many local minima and it was rather difficult to draw general conclusions on how the non-optimized property behaves with the connectivity.
We also implemented a multi-objective optimisation algorithm, in which a perturbation in ε ij was only accepted if it increased simultaneously both the explorability and the dynamical robustness of the network. However, this method worked only for tree-like networks with one or two additional links (see red dot in Fig. 4), while it was totally inefficient for more dense structures. Slightly better (but still suboptimal) results were obtained when both tasks were optimised one at a time in alternating periods.
In conclusion, sparse networks provide quasi-optimal values for both the explorability and dynamical robustness without fine-tuning many of the interaction strengths.
Self-similarity. Finally, we proved that the property of sparsity is self-similar, since on aggregating sparse interacting communities, we obtained larger sparse communities (see Supplementary Information for details). For example, joining two networks with a tree-like topology using a single link led again to a network with a tree-like topology. Similarly, if sparse networks with S nodes have − aS b links, with a and b integer constants, then joining two such networks with S and S′ nodes using b links leads again to a sparse network with + ′ − a S S b ( ) links. Therefore, the optimal features of sparsity are conserved on assembling or disassembling processes, thereby avoiding any drastic change in the stability 37 .

Discussion
Our approach provides a theoretical insight into why sparsity is an observed common feature in living interacting systems. Sparse networks generally offer optimal values of both explorability and dynamical robustness, whereas denser networks can only perform better if interactions are selectively tuned. Nevertheless, we observed that finding dense optimal networks with higher values of both explorability and dynamical robustness was barely feasible due to the multiplicity of the parameters that must be simultaneously tuned. Moreover, typically, the final networks have values of explorability and dynamical robustness comparable to those achievable for tree-like networks structures without the need to tune any parameters.
The results presented support the idea that sparsity is an emergent pattern of living interaction networks and this has implications for the understanding of the relationship between stability and complexity in real ecosystems. Indeed, sparsity may play a key role in the resolution of the so-called complexity-stability paradox 38,39 , in which highly biodiverse ecosystems will probably be unstable. The essence of the argument 38,39 , can be summarised as follows. The linearised dynamics for the population density around a stationary state depends on what is known as the community matrix, M. If all the eigenvalues of M have negative real parts, then the stationary point is also stable against small perturbations of the stationary populations. A null model corresponds to assume that The optimal tree-like network (containing S edges) is represented by the orange dot. The red dot in the right top corner represents a tree-like network with one additional link (i.e. a network with S + 1 edges), obtained using a multi-objective optimization for both the explorability and the dynamical robustness; such a method appeared to be inefficient for higher values of the connectivity. Each network was optimized using a stochastic hill climbing algorithm 35 with parameters T = 200 and δε = .
M is a random matrix with diagonal elements (the self-interactions) equal to −d < 0, whereas the off-diagonal elements are zero with probability 1 − C and with probability C are drawn from a probability distribution with zero mean and variance σ 2 . Under this null hypothesis, one finds (see 38,39 , and references therein for rigorous results) that the stationary point is unstable with probability 1 if σ > CS d, where S is the number of species in the ecosystem (a measure of its biodiversity). This result holds if S is assumed to be large enough. Thus if σ and d do not have a peculiar scaling with the network size, highly complex ecosystems (i.e. with high CS) are not stable: a prediction in contradiction to empirical data 11,37,40 , However if the interaction network is sparse, i.e. ∼ C S 1/ , the above inequality becomes independent of S and the stability of the ecosystem is not threatened by high biodiversities: sparsity in ecological interaction networks allows for stable large living interacting systems 11,19,37 , As a matter of fact, such a scaling relationship is supported by the empirical observation (see Fig. 1).
Recent theoretical findings show that an increase in the interconnectivity between multiple systems composed themselves of interacting units can have a strong impact on the vulnerability of the whole system 41 . In the same vein, our results provide a theoretical understanding of this feature. We suggest that sparsity is a key feature allowing living systems to be poised in a state that confers both robustness and adaptability (explorability) to best cope with an ever-changing environment and to promptly react to a wide range of external stimuli and to resist to perturbations.
Our ideas could be applied to understand the emergence of sparsity in some non-biological systems, as it has been empirically observed in human-built networks 42 . The underlying hypothesis of our approach is that living system interactions, unlike physical interactions (e.g. electromagnetic interactions among charged particles), can evolve/adapt to turn themselves on or off. The same idea holds for many non-biological networks, i.e. interactions can be selective and change in time, although the equations governing the underlying dynamics can be unknown.
Finally, we stress that our results do not depend on the specific details of the system and thus can be applied in many other fields. For example, a possible application might be in the design of artificial learning machines, such as deep neural networks 43,44 , There is mounting evidence that deep learning often finds solutions with good generalisation properties 43,45 , and it has been shown recently 46 that to achieve such a good performance, it is crucial to have regions of the optimisation landscape that are both robust and accessible, independent of the particular task or of the training data set. On the other hand, maximisation of computation efficiency is a crucial point when designing learning machines: deep networks are very dense as each node is connected to all other nodes of the adjacent layers 44 , which makes multilayer neural networks computationally hard to train. Our solution suggests that designing sparse neural networks will increase the explorability of the system while improving the convergence and robustness properties of the existing optimisation algorithms.