The ambiguity of nestedness under soft and hard constraints

Many real networks feature the property of nestedness, i.e. the neighbours of nodes with a few connections are hierarchically nested within the neighbours of nodes with more connections. Despite the abstract simplicity of this notion, various mathematical definitions of nestedness have been proposed, sometimes giving contrasting results. Moreover, there is an ongoing debate on the statistical significance of nestedness, since random networks where the number of connections (degree) of each node is fixed to its empirical value are typically as nested as real ones. By using only ergodic and unbiased null models, we propose a clarification that exploits the recent finding that random networks where the degrees are enforced as hard constraints (microcanonical ensembles) are thermodynamically different from random networks where the degrees are enforced as soft constraints (canonical ensembles). Indeed, alternative definitions of nestedness can be negatively correlated in the microcanonical one, while being positively correlated in the canonical one. This result disentangles distinct notions of nestedness captured by different metrics and highlights the importance of making a principled choice between hard and soft constraints in null models of ecological networks.

Scientific Reports | (2020) 10:19903 | https://doi.org/10.1038/s41598-020-76300-1 www.nature.com/scientificreports/ properties of the real system, in order to have a tailored benchmark. Null models were used, for instance, to detect the effect of the degree sequence [22][23][24][25] . Actually, to properly define a null model, the approaches to follow can be, essentially, 2: 1) to impose constraints exactly or 2) to impose constraints on average, respectively microcanonically and canonically, according to the Statistical Physics jargon. Theoretical tools from statistical physics are not new for the analysis of ecological system: they are commonly used to investigate patterns in biological networks, targeting, from time to time, hierarchical systems 26 , bipartite structures 27 and topological properties of scale-free networks 28 .
Regarding the former case, beside various approaches, the algorithm of Ref. 29 only was shown to be ergodic 30 , i.e. to visit uniformly the phase space and, thus, to provide unbiased predictions. Instead, in the canonical approach, constraints are satisfied on average and, due to its derivation, is ergodic a priori 4,31 . The canonical approach allows for some noise in the data: indeed if there is some noise, an existing pollinator-plant interaction may not be detected and the microcanonical approach will not consider the real configuration among the possible ones, while the canonical one will. Using the canonical approach, Ref. 32 compared the metric introduced in Ref. 15 with a null model preserving the degree sequence and found that in most of the cases the degree sequence is responsible for the high value of the nestedness (actually, the null model implemented in Ref. 32 is out the regime of validity). The recent contribution of Payrató-Borràs et al. 33 came to similar conclusions, using an improved null model still preserving the degree sequence, but valid for any level of link density of the network, although using an approximated formula for the average of NODF in the ensembles; subsequently the same group enlarged their analysis to a wider number of nestedness metrics 34 .
As presented above, in the literature, several papers compared the measurements with various null models [22][23][24][25][32][33][34][35] , but rarely the issue of ergodicity was targeted. For the first time here, we investigate the differences of the micro-and canonical approach in discounting the degree sequence for the analysis of the nestedness. Let us remark that both the null models implemented are ergodic, i.e. they explore the phase space homogeneously, in order to have unbiased and unequivocal results.
In the present paper, we provide an example in which two null models ergodically discounting the same information, i.e. the degree sequence, display opposite correlations between two among the most used nestedness measures. Beside providing another example of the statistical ensemble inequivalence, the main message of our manuscript is that choosing to quantify the amount of nestedness is a subtle task that has to be carried out carefully. Indeed, only being aware of the behaviour and the peculiar properties of the various approaches and options permits to derive the right conclusions from the analyses: our paper provides the necessary knowledge to handle properly the study of the nestedness of a real system. In this sense, we do not provide any univocal indication on which nestedness definition measure should be used or on which is the correct way to discount the information encoded in the degree sequence: both the definitions analysed have their own sense and both the null models examined satisfy their own rationale. Nevertheless, it is crucial to know the properties of the various tools that one is handling in order to derive the proper conclusions from the nestedness analysis of a real system.

Methods
A bipartite network is defined by two sets of nodes L (of size N L ) and Ŵ (of size N Ŵ ) called layers and by the prescription that connections are allowed only between the layers and not inside them. Thus, a bipartite network can be univocally described by its biadjacency matrix M , i.e. an (N L × N Ŵ )-matrix, whose entries m iα = 1 if a link exists between i ∈ N L and α ∈ Ŵ and m iα = 0 otherwise. We will call a network perfectly nested (PNN, Perfectly Nested Network) if for every pair of nodes i, j belonging to the same layer with degrees d i , d j , if d i ≤ d j then all neighbours of i are also neighbours of j. This type of network is also called chain graph 20 or double nested graph 36 . In the following we shall use the previous definitions for generic biadjacency matrices and related quantities, but we shall add an asterisk * whenever considering quantities measured on real networks. Nestedness measures. NODF. One of the most popular measure of nestedness, namely the Nestedness as a measure of Overlap and Decreasing Fill (NODF) was introduced in 2008 by Almeida-Neto et al. 19 . Such measure is based on the overlap between the neighbourhoods of nodes with different degrees. Given a generic bipartite graph G Bi , the NODF expression reads where K = [N L (N L − 1) + N Ŵ (N Ŵ − 1)]/2 is a normalisation factor to let the measure go from 0 to 1, k i and h α are respectively the degrees of node i ∈ L and α ∈ Ŵ , and θ is the Heaviside step function with the convention θ(0) = 0 . The step function ensures that the overlap is only counted when the degrees of the nodes are different and that the denominator is the minimum of the two vertices' degrees.
Stable-NODF. Due to the instability of the previous measure, with respect to small fluctuations on the degrees of the nodes, another version was proposed in 5 . The difference relies in considering also the contributions coming from couples of nodes with equal degrees; we will call it stable-NODF or sNODF. It is calculated as www.nature.com/scientificreports/ where K is the same normalisation factor as in Eq. (1) and the denominator this time is the minimum between the two degrees, that in Eq. (1) was guaranteed by the theta step function.

Spectral nestedness.
A recently proposed measure of nestedness 21 considers the spectral radius of the network, i.e. the largest eigenvalue of the adjacency matrix. We will thus call it spectral nestedness (SNES). The adjacency matrix A of a bipartite network can be expressed in terms of the biadjacency matrix as where M T is the transpose of the biadjacency matrix M and 0 N×N is a N × N-matrix whose elements are all zeros. Note that the adjacency matrix of the network is symmetric, yielding all real eigenvalues. The definition is based on two main theoretical results: • The bipartite network that has the maximum eigenvalue in the set of connected networks with given n nodes and L links is a perfectly nested network 36 ; • Among all bipartite networks with a given degree sequence on one of the two layers, the one that maximises the spectral radius is the PNN, defined at the beginning of the present section 20 .
Normalised spectral nestedness. The spectral radius, though, has a strong dependence on the size of the network and on its density. It is well known that the maximum eigenvalue of a bipartite network with L links is bounded from above by √ L and that the only network for which (M) = √ L(M) (if it exists) is a complete bipartite network 20,36 .
For this reason we decide to introduce nSNES, where we normalise the measure with the square root of the number of edges: Although the nSNES ranges from 0 to 1, the drawback of this normalisation is that a perfectly nested matrix that is not full will not have a perfect score of 1.
Null models. In the present paper, we aim at understanding the role of the degree sequence in the formation of bipartite nested structures. Thus, we would need a sort of network benchmark with the same degree sequence, but otherwise maximally random. This approach has strong similarities with Statistical Mechanics: actually, the recipe is to build an ensemble and fix the node degrees on it. As in the standard Statistical Mechanics, those constraints can be imposed on average, as in the canonical construction 4,31,37-39 , or considering stricter constraints, as in the microcanonical formulation 29,30 . The two approaches are known to be non equivalent [40][41][42][43][44][45][46] and indeed such non equivalence is going to be crucial in the following.
After the randomisation with the null models, our aim is to quantify the statistical significance of the measures by computing the z-scores of the measures, calculated as where σ X is the standard deviation and X the considered quantity.
The canonical approach: the bipartite configuration model. The Bipartite Configuration Model (BiCM 47 ) is the bipartite extension of the entropy based null model 4,31,37-39 . The strategy is inspired by work by Jaynes 48 , which derived the canonical ensemble of Statistical Mechanics from Information Theory principles. The recipe is pretty simple: first, define an ensemble of all possible physical configurations, and then maximise its Shannon entropy constraining the relevant information about the system (in the case of Information Theory, the energy): the result is exactly the canonical ensemble. The maximisation of the Shannon entropy represents the crucial step: it can be interpreted as assuming maximal ignorance about the the non constrained degrees of freedom of the system.
Following the same strategy, starting from a real network, we can define M the ensemble of all possible biadjacency matrices with the same number of nodes (nodes represent the volume in Statistical Mechanics). The Shannon entropy associated to the ensemble is S = − M∈M P(M) ln P(M) and we can maximise it, constraining the degree sequence (Note that while in the Jaynes derivation of the Statistical Mechanics, the constraints, i.e. the energy, was a global one, the degree sequence represents a local one. Actually, the local constraint is responsible for the nonequivalence of the microcanonical and canonical ensembles.). The entropy maximisation leads to an exponential probability for a generic biadjacency matrix M: where C(M) is the vector of constraints and θ the associated Lagrangian multipliers 37 . At this level the formula 4 is just formal, in the sense that the value of the Lagrangian multipliers is unknown. At the end of the day, we want a "tailored" benchmark for our real network, i.e. something with the same degree sequence, but otherwise www.nature.com/scientificreports/ completely random. In this sense, it is natural to maximise the likelihood of the real network in order to get the value of θ 38,39 . If � C(M * ) is the value of C measured on the real network, the previous condition is equivalent to . The exact solution for the probability P(M) can be factorised as the product of probabilities per possible link: where p iα is the probability of existence of the link connecting nodes i and α . Let us remark that the factorisation (5) is possible only when the constraints are linear in the biadjacency matrix. For other nonlinear contraints, the probability per link may not be analytical and other methods are necessary to obtain the probability per graph (see for instance 49 ).
In the case of the BiCM, p iα is a function of x i and y α , which are simple reparametrisations of the Lagrange multipliers associated to the observed degrees ( k i and h α respectively): Their numerical value is determined by solving the likelihood-maximisation equations: k i and h α being the degree of the node i and α respectively.
The microcanonical approach: the Curveball algorithm. The microcanonical approach, differently from the BiCM, keeps the degrees of all nodes in the system constant. In a sense, it has a stricter ensemble (just all configurations with the given degree sequence are allowed) and all allowed configurations have the same probability. Such approach is computationally costly since the probabilities of links in the system are not pairwise independent and the fastest way of spanning the ensemble of networks with a given degree sequence relies on swapping endpoints of links iteratively. In the present manuscript, the ensemble was sampled using the strategy of 29 .
We will refer to this model as Curveball, as in the original paper; in 30 it was shown that such approach is ergodic.

Results
In this section we are going to present the results of our analyses on artificial and real networks. To test the measures and models, we analyze a set of 40 pollination networks taken from the Web of Life dataset (www.web-of-life. es). They represent ecological mutualistic networks of plant-pollinators. From the set of all available pollination networks in Web of Life, we selected only the binary ones, in order to avoid issues regarding binarisation. All of the considered networks are generally of small size, the smallest being of only 20 nodes while the biggest one consists of 1500 nodes. The density of the networks varies between 0.01 and 0.5. For the sake of completeness we remark that only 24 out of 40 networks of our dataset are actually made of a single connected component, the other including few disconnected components with more than one node. In the following, we compare the various measures and state their significance respect to the various null models.
Measure differences. First, in order to study the behaviours of the previous measures, we compare them on the above-mentioned dataset. Fig. 1 shows that indeed the normalised SNES is highly correlated with NODF (actually, it is not true for the non normalised version of the spectral nestedness, due to its dependence on the total number of link). In a sense we may think that indeed, while they differ in the philosophy, the two measures are capturing the same structure, as stated in 21 . After a detailed comparison with the appropriate null models, we will see that it is not the case.
Degree sequence versus nestedness. The degree sequence of the network carries some information about the nestedness of the system, the extreme case being the Perfectly Nested Network (PNN in the following) one. Actually, in this case, the degree sequence identifies completely the network and both the micro-and the canonical ensembles are composed by a single network, i.e. the PNN one. This was already observed in 50 for the microcanonical ensemble, but it is surprisingly true also for the canonical ensemble. While the technical details can be found in Section 2 in the Supplementary Information, a pictorial representation is shown in Fig. 2.
Thus in these cases the degree sequence captures the level of nestedness of the whole system, and the statistical significance of the measure loses any value. Actually, even when the network is close to a perfectly nested one, its configuration model ensembles contain a limited number of configurations, that furthermore are all highly nested networks. Thus, a real network may show a high value of the nestedness measure (whatever it is), which is, nevertheless, statistically non significant with respect to a null model discounting the degree sequence: actually in such a case the high value of the nestedness is already captured by the degree sequence. We will examine in more details the role of the null model in the following sections. Actually, there is another issue generating a bias the canonical ensemble. Indeed, a network sampled from the ensemble can present some isolated nodes that do not contribute to the measurements of both NODF and SNES (and their modifications). Given the steep power law degree distribution of many of the considered networks, this will typically be the case. For further details, in Section 3 of the Supplementary Information we analyse the frequency of isolated nodes in generating the canonical ensemble. We will discuss this issue and how it generates a bias in greater detail focusing on each of the two measures in the following paragraphs.
NODF versus null models. The displacement of the NODF measures between the two ensembles is the result of multiple effects. The most evident bias is caused by the normalisation factor that is the denominator in Eqs.
(1) and (2). A network sampled from the BiCM ensemble will have, on average, the same number of links of the original network, but many isolated nodes (see the Supplementary Information Section 3 for more details). This is due to the small link probabilities related to nodes of low degree in large networks, that sometimes give rise to an empty row or column in a sampled matrix. These nodes, therefore, do not contribute to the total NODF or sNODF, and one has to choose how to handle the normalisation factor in 1 and 2.
If one chooses to consider the number of connected nodes of the sampled network (as in the original definition of the NODF), this will generate a positive bias by having a comparable quantity divided by a lower denominator (the number of the connected nodes in each realisation can be only smaller than the value of the real network). We call this approach heterogeneous.
Otherwise, considering the normalisation factor of the original network will introduce contributions even from isolated nodes, thus altering the philosophy of the original definition. Moreover, such approach will introduce a bias in the opposite direction, dividing by a factor that is larger than what it should be if considering only the connected network. We call this approach homogeneous.
Both choices are equally admissible, depending on the interpretation of the comparisons one wants to follow. Personally, we think that the normalisation should not involve the isolated nodes, as in the original definition, i.e. we prefer the heterogeneous normalisation. For completeness, in the next subsections we will consider both of them. Interestingly enough their differences do not affect the conclusions.
On top of this, another effect to be considered is the presence of fluctuations in considering the degree sequence. As mentioned in the previous section, both ensembles contain only one configuration in the case of a perfectly nested degree sequence and their measures are trivially exactly the same. When the two ensembles separate for a non-perfectly nested matrix, the canonical ensemble produces some variance in the degrees of the nodes. This effect is not present in the microcanonical ensemble, where the degrees of all nodes are fixed deterministically. Such an effect has an impact on the NODF and in particular it provides new evidences regarding the non equivalence of the various ensembles.
SNES versus null models. We observe that the spectral radius is slightly overestimated in the canonical model. Our guess for this behaviour is that on average, out of two matrices with the same number of links, the one with the smallest number of nodes has the largest radius, so when a sample of the canonical model has an empty row or column it has, on average, a higher radius. Some evidences for this behaviour are given in the Supplementary Information Section 4. Still, we are not able to evaluate such discrepancy. Regarding the nSNES, the overestima- www.nature.com/scientificreports/ tion seems to be strongly increased by the normalisation factor. Although this is intuitive, since � √ L� < √ �L� because of the non-zero variance of L, the fact that the spectral radius has an intrinsic dependence on L makes it hard to evaluate precisely.

The significance of the nestedness measures with respect to the various statistical ensembles.
Bearing in mind all of the considerations of the previous paragraphs, we can interpret the z-scores of the table in Fig. 3. The four canonical NODF columns of the table refer to the z-scores of the two variants of NODF, with the two different normalisations with respect to the canonical ensemble. A representation of the differences among the models and measures is also given in Fig. 4.
In the case of the heterogeneous normalisation the z-scores are all negative, because of the overestimation of both NODF and sNODF in the ensemble. There are, though, important differences between the NODF and sNODF measures in some cases, which are mainly due to the presence of many nodes with the same degree.
For the homogeneous normalisation, there is still a certain agreement in the signs of the z-scores of NODF and sNODF, with the same caveat discussed above for the heterogeneous case. In opposition to the heterogeneous columns, the z-scores are positive in most of the networks analysed, in agreement with the discussion of "NODF versus null models" section). In this sense, it is striking that the choice of the normalisation factor may drive to opposite conclusions, regarding the statistical significance of the measure on the real network.
Then we have the columns of SNES: similarly to the heterogeneous normalised NODF, even in this case, the canonical null model has all negative z-scores, due to the slight overestimation of the SNES. The second-last column contains the microcanonical SNES z-scores. Figure 2. An example of a perfectly nested network with its probabilities per link from the BiCM: at the first step, the first row and column are full, and the degree is respectively 12 and 8. So the link probabilities must be exactly one, for preserving the row sum and the column sum. At the second step, since the last row and column have degree 1, the remaining entries must sum to 0, yielding all zeros. Again, at the third and fourth steps the rows and columns that are completely full or empty univocally determine the respective probabilities to be 1 or 0. At the end of this process, the link probabilities are all set to 0 or 1, so the corresponding canonical ensemble contains only one matrix.
Scientific Reports | (2020) 10:19903 | https://doi.org/10.1038/s41598-020-76300-1 www.nature.com/scientificreports/ As it can be observed from the matrix, there is no agreement between the column of the SNES and the sNODF columns in the microcanonical ensemble. Note that for the microcanonical null model, the z-scores of the nSNES were not reported since they correspond to the one of the SNES (the normalising factor cancels). A hint is given by the assortativity z-scores, whose Pearson correlation coefficient with the SNES scores is 0.84, while SNES and NODF anti-correlate with a score of − 0.88.
In order to investigate this difference, we generate a scatter plot of the realisations of the different ensembles (Fig. 5), plotting the NODF against the SNES of the sampled networks. The results are striking: the two measures are highly anti-correlated on the microcanonical ensemble, while this effect is hindered by the fluctuations in the canonical ensemble.
Using the other proposed measures, the results are always similar when comparing a NODF measure and a spectral nestedness measure. NODF and SNES are actually capturing different ways of being "nested". This is easily seen on a synthetic very small network, of size 8 × 9 . We generate a sample from the microcanonical model and see how the matrices maximising NODF and SNES are made (Fig. 6).
Actually the NODF-maximising matrix has one of the smallest value of assortativity, while the one that maximises the SNES presents a big hub of the highest degree nodes and two smaller disconnected subgraphs, see Fig. 6. Roughly speaking, on the one hand, the SNES prefers networks in which highly connected nodes link to highly connected nodes, since they are sort of carrying the "mass" of the adjacency matrix (which is what the www.nature.com/scientificreports/ spectral radius is measuring). On the other hand the NODF, due to the denominator of its contributions, prefers to link poorly connected nodes with highly connected ones, thus focusing on disassortative configurations. Regarding the anti-correlation between the NODF (or similar definition) and the assortativity, other studies got to similar conclusions 32,52 ; as far as we know, there were no evidences regarding the opposite behaviour of the SNES. Let us underline that the (anti)correlation between the NODF or SNES and the assortativity is present only when discounting microcanonically the degree sequence: in real data, such correlation is not evident, as Fig. 7 shows. Otherwise stated, the (anti)correlation is present only when the contribution of the degree sequence is discounted. More details regarding the correlation between assortativity and the various nestedness metrics can be found in the Supplementary

Discussion
While the abstract idea of nestedness in networks is quite straightforward, we saw that a mathematical definition capturing the degree of nestedness of a real system is much less trivial. As a consequence, while nested structures are ubiquitously observed across several networks, measuring the actual level of nestedness along with its statistical significance remains a challenging task.
In the present manuscript we investigated in details different metrics of nestedness in both real-world and synthetic networks. In particular, we mainly focused on two measures, NODF 19 and SNES 21 , and some of their modifications 5 . When applied to real networks, these metrics go in the same direction, as they give positively correlated results. www.nature.com/scientificreports/ We then moved to discount the contribution of the degree sequence to the different nestedness measures. Literally, according to the case of study and to the available information on our system we can create suitably chosen series of randomised copies of our graph (ensembles). This procedure allows us to use the machinery of Statistical Physics to assess the significance of our measurements.
Thus, for our aim, we can define null models preserving the degrees of nodes either as hard (microcanonical ensembles 29 ) or as soft (canonical ensembles 4,31 ) constraints. Otherwise stated, we are using the extensions of the microcanonical and canonical ensembles to complex networks in order to discount the information carried by the node degrees: the degree sequence is supposed to have an effect on the nestedness 32,33 , thus we want to focus on the information carried by the different metrics that cannot be explained by the degree sequence only. Let us remark that the null models implemented are ergodic, i.e. they explore the phase space uniformly.
First, we concentrated our attention on Perfectly Nested Networks (PNN). A PNN has a degree sequence that admits only a single network, i.e. the PNN itself, irrespective of whether the degrees are treated as hard or soft constraints: both the microcanonical and canonical ensembles of a PNN are composed by the PNN network only. Otherwise stated, there exist perfectly nested degree sequences and each of them defines univocally a single network, i.e. the PNN one. In the case of PNNs, thus, the value of the nestedness is completely due to the degree sequence only. But what happens when the network is not perfectly nested?
We compared the values of NODF and SNES measured on real networks with the expectations of, respectively, the microcanonical and canonical ensembles. As theoretically demonstrated in other studies [40][41][42][43][44][45][46] , the two ensembles are not equivalent, thus they should be characterised by different macroscopic properties. Literally, we found that the two families of nestedness metrics are negatively correlated when the microcanonical ensemble is used, while they are positively correlated in the canonical ensemble. Actually, the fluctuations of the canonical ensemble cover the real behaviours of NODF and SNES. Instead, once the degree sequence is fixed as a hard constraint, the level of nestedness is influenced by higher-order correlations between the degrees themselves, and in particular the assortativity of the network. Indeed, the two classes of measures of nestedness give different results in the microcanonical ensemble, when considering networks with different assortativity: NODF tends to give larger values of nestedness when the network is disassortative, while SNES tends to give larger values of nestedness when the network is assortative. Otherwise stated, we present an example in which the same information, i.e. the degree sequence, is ergodically discounted, but the sign of the measured correlations are opposite if compared with a microcanonical of with a canonical null model.
Thus, other than the choice of the measure, if checking for the statistical significance of the nestedness of a system, one should make a principled choice of the ensemble used as a null model in the analysis 45,46 . The microcanonical ensemble, which treats degrees as hard constraints, should be preferred if the observed degrees are error-free, i.e. if they are the actual values of the property to be kept fixed in the null hypothesis. If one suspects that the observed degrees are instead subject to some sort of error (e.g. measurement errors, incomplete data collection, poor sampling, etc.), then the microcanonical ensemble should be avoided, as it will give zero probability to the true (undistorted) configuration and to any configuration with the same degree sequence as the www.nature.com/scientificreports/ true configuration. In this case, we suggest to use the measures that present the smallest biases for fluctuating degrees, i.e. the SNES and the homogeneous sNODF or NODF. Let us finally remark that our paper does not provide any indication on which is the nestedness metric that should be used, or on which is the right null model to be implemented in order to state the statistical significance of the nestedness measured. In a sense, each nestedness measure has contraindications, and every null model, even if discounting the same information, has its peculiar properties. In this sense, it is crucial to know exactly the behaviour of the ingredients we are handling. Our contribution is in highlighting odd behaviours, previously undetected, that, if not under control, can take to unjustified conclusions.
Nevertheless, in light of our results, the question remains on how to tackle the problem of the nestedness in real system, even after a proper and justified choice of the nestedness measure. An easy solution can be, once the null model has been chosen, to report for the chosen nestedness measure, both the average over the ensemble and the z-scores on the real network. The former value provides an evaluation of the nestedness as encoded by the degree sequence, the latter how significant is the observed nestedness, once the degree sequence is discounted. Figure 6. Top: a sample from the microcanonical configuration model ensemble with the relative scores of nSNES and sNODF. Bottom: the same sample with scores of sNODF against assortativity (left), nSNES against assortativity (right). Different colours are used for sampled networks that result connected or disconnected. The highlighted networks have high values of nSNES or sNODF, and show that for their extreme values, the systems can be disconnected (left) or barely connected (right). The left matrix, with a high nSNES, has a really assortative configuration, while the right one has a high sNODF and is highly disassortative. We do not exclude disconnected networks in our analysis since it could be a possible configuration for an ecological system.