Consider a system that evolves as (1). We state our main result in the simplest formulation of a directed unweighted network M, which means that M ij = 1, if the change of i depends on the change of j. In such systems, the normalized vector \({x}_{i}={X}_{i}/\sum _{j}\,{X}_{j}\) “quantizes” if there is one single cycle remaining in M. This means that all entries in x i are multiples of some minimum value x min . This is captured in the following theorem:

Theorem 1 (Eigenvector Quantization). Let M be a binary matrix with entries M ij ∈ {0, 1} and diagonal entries M ii = 0 for all i ∈ {1, …, N}. Let G be the directed network with directed adjacency matrix M. Let X(t) = (X 1 (t), …, X N (t)) be an N-dimensional state vector, whose components X i (t) evolve according to (1). For all initial conditions X(0), except for a set of Lebesgue-measure zero, the normalised vector x(t), defined component-wise as \({x}_{i}(t)={X}_{i}(t)/\sum _{j}\,{X}_{j}(t)\), converges to a stable fixed point x:= lim t→∞ x(t), for which the following holds:

Eigenvector Quantization: Suppose G contains only one single cycle. Then any component x i can be expressed as

$${x}_{i}={n}_{i}{x}_{{\rm{\min }}},$$ (2)

where x min is the minimal non-zero component and n i is a natural number. The value of x min is taken by the cycle-nodes, and the integer n i ≥ 0 is the number of directed paths that lead from cycle-nodes to node i. If there are no paths from cycle nodes, then x i = 0.

For the proof of the theorem, see (SI). Note that the theorem gives a sufficient condition for the eigenvector quantization of the node states. A graphical illustration of the theorem is shown in Fig. 1. We compare two networks that are identical except for one link. The network in Fig. 1 on the left-hand side contains two cycles, indicated by blue arrows, the network in Fig. 1 right-hand side has only one last remaining cycle, again indicated by blue arrows. The node color corresponds to the respective components of the normalized state vector, x, rescaled by the minimum value, x min . The histogram to the right of each network figure indicates the number of nodes, whose components share a given value of x/x min . It can be seen from the histogram in Fig. 1 right that the state vector takes only discrete, equidistant values when the network has one remaining cycle. Once the quantization occurs, one can immediately determine the number of paths by which a node can be reached from the cycle, x i /x min . Node A can be reached via two paths from the cycle, while node B can be reached via four. For the multi-cycle case (a), the number of paths no-longer coincides with the states.

Figure 1 Graphical demonstration of the quantization theorem. Directed networks M containing two cycles (left) and one cycle (right). Cycles are in the shaded area. The node color indicates the state, x i , (value of the component of the state vector, x) in units of the minimal value x min . The histograms show the number of nodes in a given state. In the right-hand figure we see the quantization of states due to the presence of a single (!) remaining cycle. The occurrence of quantization at the last remaining cycle can be used as a precursor signal. The state x i /x min in the single-cycle network (right) coincides with the number of directed paths from the cycle to node i. Node A can be reached through two paths from the cycle, while node B can be reached by four. For the multi-cycle case (left), the number of paths no-longer coincides with the states. Full size image

The key mechanism behind quantization can be explained as follows. The state of a given node does not only depend on the nodes which feed directly into it, but through them, indirectly also depends on a large number of other nodes, and possibly the entire network. This global interdependence is best understood through the equation that is satisfied by the stable fixed point of the normalized state-vector x = lim t→∞ x(t), namely the eigenvalue equation (Lemma 1 in SI):

$${\lambda }_{1}{x}_{i}=\sum _{j}\,{M}_{ij}{x}_{j},$$ (3)

where λ 1 is the maximal eigenvalue of M. This equation couples all node states. It is a well-known fact that λ 1 = 1 if the network contains only one cycle34. In this case \({x}_{i}=\sum _{j}\,{M}_{ij}{x}_{j}\). The state of any node without incoming links must therefore be zero, and any node-state that receives only in-links from such zero-nodes also vanishes. Consider, on the other hand a single cycle in isolation. Each of its nodes inherits the states of their respective parent nodes, and therefore all cycle-nodes must have the same non-zero state. Now, if the cycle is embedded in a network, it has outgoing and incoming links. The incoming links cannot originate from the cycle, or else there would be nodes that are part of more than one cycle. Hence, they must come from vanishing nodes and cannot contribute to the cycle states. Therefore, also the nodes in the non-isolated cycle still have the same states. Finally, consider a node that has multiple in-links. Its state is simply the sum of the node-states pointing to it. Consequently, if a node has two, three, or four paths leading to it from the cycle (see node A or B in Fig. 1 right), its state is two times, three times, or four times the state of a cycle-node.

The magnitude of state x i also captures the eigenvector centrality of the corresponding node i in the network. The eigenvector centrality is either defined iteratively or directly via Eq. (3), with λ 1 equal to the largest eigenvalue44. The centrality x i does not just measure the number of neighbors that i receives links from, but also accounts for their importance. So, a node with a high eigenvector centrality is one that receives links from very central nodes. The definition depends on the convention of the arrows for a given adjacency matrix. Equation (1) therefore gives rise to the eigenvector centrality in a dynamical way. Likewise, one may ask to which extent the quantization of the state vector changes for other centrality measures.

In practice, it is more useful to consider a reverse direction of the theorem, i.e., to predict one cycle from the eigenvector quantization. Indeed, the reverse direction is not true in general. When λ 1 = 2, 3, …, it may happen that we observe the quantization too. However, for sparse networks, where λ 1 is typically close to one, the reverse is true. In particular, this is the case for the Jain–Krishna model, which is presented in the following. In this case the eigenvector quantization can be used to determine the structural stability of and can provide a warning signal for impending collapse. Further, one may determine the path-distance of a given node to the last remaining cycle. In the following, we discuss several examples in which the quantization can be detected. In the first example, we focus on the aforementioned Jain–Krishna model where we demonstrate explicitly how well it works to actually predict a collapse.

Co-evolutionary complex systems — Jain–Krishna model

We consider a simple model for the evolutionary dynamics of a toy ecosystem33. It was used to explain the rapid changes (punctuated equilibria) in species diversity and the network mechanisms behind it. In particular, the model shows large-scale collapses in diversity. We briefly describe the model, explain how our main result applies, and show how predictions about the average time to collapse can be obtained.

The Jain–Krishna model is a co-evolving network model that distinguishes between a fast time scale, at which populations evolve on a fixed interaction network M, and a slow time scale, at which the network evolves through selection and mutation of its species. The population X i of species i is a proxy for its abundance. It evolves according to the linear catalytic dynamics given by (1). The relative abundance x i (t) converges to a fixed point x i . Thus, x i is endogenously defined by the interaction network M and can be used as a measure for the fitness of species. The network is updated on the slow timescale as follows: the species with the smallest value of x i (along with all its links) is eliminated. A new species is introduced, again labelled by i, which forms relations (links) with, on average, m randomly selected already existing species j. This is realised by independently assigning in-links to (M ij = 1), and out-links from (M ji = 1) the new species, i, both with the same probability m/(N − 1). In33, a range of m from 0.05 to 2 is considered. For m < 1, the model shows collapses. The Jain–Krishna model depends on two characteristic parameters: the average connectivity m, and the maximal diversity N. The observables in the model are the abundances x i and the diversity given by the number of species that have a positive abundance,

$$S=\mathop{\sum }\limits_{i=1}^{N}\,\Theta ({x}_{i}),$$

where Θ is the Heaviside step function that equals 1 whenever its argument is greater than 0.

We briefly describe the mechanisms that lead to species diversification and to collapse. If there are no cycles in the network, most species must die out. We show this formally in the SI (case λ = 0 in Lemma 1), but here we provide an illustrative example. Consider only two species connected by a directed link. The receiving species grows proportional to the giving species. Asymptotically the former species gains a positive relative abundance, whereas the latter vanishes in relative terms. This argument can be extended to all cycle-free graphs, for which it can be shown that the only the species at the end of chains attain a positive relative abundance and the rest have asymptotically vanishing relative abundance (Lemma 1 of the SI). If there is at least one cycle, many species have positive abundance, they exist. Existing species are either part of a cycle, or have incoming paths from a cycle. In the network updating process the system may change, if a hitherto unfit species dies out and upon its re-introduction attaches itself to a fit (very abundant) species, thus becoming fit itself. The network updating process can also lead to the creation of new cycles, which can lead to a strengthening of the existing system33. Eventually, cycles can also break, whenever cycle-nodes become the least abundant species. This inevitably happens, since those species which are not part of a cycle either are already well adapted by influx from highly abundant species, or keep being replaced until they become well adapted. When all cycles break, all self-sustaining structures disappear, and the ecosystem collapses.

Every collapse is always preceded by a phase where only one remaining cycle exists in the network, see (SI). This is demonstrated in Fig. 2 left, where a typical example of the Jain-Krishna collapse mechanism is shown. The detection of this last cycle thus becomes a warning signal for the impending collapse. Since the population dynamics is governed by (1), the quantization theorem applies for the Jain–Krishna model. Whenever the species vector x quantizes a warning signal is produced; the network now has only one cycle left and is potentially very close to collapse.

Figure 2 Left: sample run of the Jain-Krishna model. ρ is the fraction of populated nodes, λ 1 is the largest eigenvalue. The time point t cycl denotes the moment of entering the one-cycle phase (where the eigenvector quantization appears). We see that the collapse is preceded by the one-cycle phase. For comparison, we add two other precursor measures — volatility spectral radius λ V and cross-correlation spectral radius λ C (their definitions are in SI). We observe no signal before the collapse, the signal appears after the collapse. Right: comparison between numerical simulations (symbols) and the theoretical prediction of (4) (dashed line) of the expected time-to-collapse, 〈T〉, for various system sizes N ∈ {25, 100, 200}, and connectivities, m between 0.05 and 0.375. Simulation results follow the theoretical prediction closely, independent of system size and m. Full size image

We can estimate the expected time-to-collapse once quantization is observed. This calculation is shown in detail in the (SI). Here we sketch the main idea: First, the probability of removing a cycle-node in the critical regime is p = p w p c , where p w is the probability that the cycle-node belongs to the set of the weakest nodes, and p c is the probability of choosing it from this set. The former can be approximated for sufficiently large sparse graphs (m < 1 and N ≫ 1) by the average network connectivity m, p w ~ m. This is true because a new species only survives in the ecosystem if it receives at least one in-link, which happens with probability 1 − (1 − m/(N − 1))N−1 ~ m. Otherwise, it remains unpopulated and becomes the weakest species with vanishing abundance and will be replaced again. The latter, p c can be expressed as the probability that a cycle of length L c is formed among the L weakest nodes. This probability assumes its maximum and hence most likely value p c ~ 1/e at L c ~ L/e, as shown in the (SI) (e is Euler’s constant).

In summary, we obtain that the probability of removing a cycle-node is p = m/e. It is now possible to estimate the average life-time of the ecosystem once it enters the critical phase. The probability that the system collapses within the next T steps after the appearance of the quantization-based warning signal is given by P(T) = (1 − p)T−1p. The expected time-to-collapse in the critical phase can then be expressed as

$$\langle T\rangle =\mathop{\sum }\limits_{T=1}^{\infty }\,T\,P(T)=\frac{1}{p}\equiv \frac{e}{m}.$$ (4)

We explain the intuition behind this result. Since m is the rate at which a populated species is chosen for extinction, the factor 1/m gives the average time that it takes to remove a populated node. Further, since there is a chance 1/e that this node is a cycle-node, one needs to repeat the network update on average e times to pick out one cycle-node. Figure 2 right shows the comparison of the average times-to-collapse 〈T〉 for numerical simulations of the model and the theoretical predictions from (4). We run the model on R for a range of connectivities m and system sizes N, always for 107 time steps. At each collapse, we record the time between the first observation of the quantization and the collapse. We average these times for all collapses in the run. The simulation results follow the theoretical predictions very well, irrespective of system size and m.

Finally, we compare the eigenvector quantization with other types of precursor signals. Well known precursor signals based on critical slowing down near the collapse point are given by the spectral radius of the lagged cross-correlation matrix λ c and the spectral radius of volatility λ V 40. For their definitions, see SI. We show the various precursor signals in Fig. 2 left, the cross-correlation spectral radius (green), and volatility spectral radius (grey). We observe no characteristic signals for the spectral radius before the collapse. Signals appear only after the collapse. We have run several independent simulation runs and never observed any signal for any of these two measures.

Epidemic spreading in the SIS model

Let us now discuss a model where traces of the eigenvector quantization can be observed. We consider the classical SIS model on a sparse, directed (static) network M. It describes the spreading of a disease without immunity. Nodes can be either infected (I) or susceptible (S). An infected node can recover at a rate, r, or transmit the infection along a directed link to a susceptible neighbour at a rate, β. According to the individual-based mean-field approximation, the time evolution of the probability, p i , that a node i is infected is given by23

$$\frac{d}{dt}{p}_{i}=\beta \sum _{j}\,{M}_{ij}{p}_{j}(1-{p}_{i})-r{p}_{i}.$$ (5)

When the average disease prevalence, p i , is small, i.e., 1 − p i ≈ 1, the dynamics of (5) becomes effectively linear and can be cast into the form of (1) by rescaling time t′ = βt, and setting Φ = r/β. Therefore, in this regime the theorem applies. However, because of these two approximations, the linearity condition is not exactly met. We expect a deviation of p from the exact eigenvector quantization, which appears in the limiting vector x of the linear dynamics (1), whenever the network contains just one single cycle.

We run numerical simulations of the SIS model to compare the situation of single-cycle with the multi-cycle networks. The probability that a given node is infected, p i , is approximated by the average time it spends in the infected state, as a fraction of the total simulation time. We draw 103 single- and multi-cycle graphs from the ensemble of Erdős–Rényi graphs of size N = 200 at the critical connectivity parameter, where the giant component would appear almost surely in the N → ∞ limit. We set β = 1, which corresponds to fixing a time-scale. Since for any finite value of r the infection dies out at some point (In the language of probability theory, this event happens almost surely, i.e., with probability one), we have re-initialized the model randomly, whenever this happened. For each network the duration of the simulation, including re-initializations, was 104 × (r + β), which ensures the same average number of events per simulation for a fixed network size. The average infection time, p i , is then recorded for every node and rescaled by a factor p min , that shifts the location of the first maximum to 1.

In Fig. 3 left we show the comparison between the rescaled values p/p min , obtained from numerical simulations, and the rescaled eigenvector x/x min for one randomly sampled Erdős–Rényi graph with a single cycle. The level structure of the infection duration is clearly visible and strongly correlated with the exact level structure of the eigenvector. In Fig. 3 right we show the histogram of p/p min generated by sampling 103 networks from the single-cycle ensemble and we compare it to the histogram of the multi-cycle ensemble in the inset. We choose a recovery rate r = 0.5. For r larger than β = 1, the infection has no chance to propagate through the cycle and the effect becomes less pronounced when compared to the r < β case, yet it is still visible. We color-coded the contributions according to the number of paths that lead to the contributing node from the cycle. One can clearly see the signature of the eigenvector quantization with four peaks in the single-cycle histogram, which indicates the average time of infection for nodes with zero, one, two, or three paths from the cycle. Quantization is not exact, as expected. Higher levels are not at integer values, but are shifted toward lower values. The reason for this is that the linearity assumption (5) breaks down. A node with two or three incoming paths from the cycle has a higher average duration of infection, because it can receive the infection in multiple ways. It cannot, however, receive the infection when it is already infected, which is reflected in the non-linear term p j (1 − p i ) in (5) that forces the value of the infection probability p i to be less than the multiplicity of paths.

Figure 3 Comparison between the rescaled state vector, x/x min , and the rescaled average time of node infection, p/p min , which is taken as a proxy for the infection probability of a given node. Left: situation for a single-cycle network with N = 200 species that was randomly sampled from the Erdős–Rényi ensemble. Circles represent the average time of infection of the respective nodes, squares are the components of the rescaled eigenvector. The strong correlation between the two is apparent. Right: histogram of the average infection times for 103 different networks of size N = 200 with single (main) and multiple cycles (inset). Note the presence of four peaks that correspond to nodes with zero, one, two, and three paths from the cycle. They are colored with cyan, blue, red, and green, respectively. No clustering around integer values is present for multi-cycle networks (inset). All networks are drawn from the Erdős–Rényi ensemble with at least one node being linked to the cycle via three paths. Full size image

Generalizations of eigenvector quantization

Here we demonstrate that the eigenvector quantization is still approximately observable for mild generalizations of eigenvector centrality of binary networks. We focus on three examples — first is the Katz centrality, second the eigenvector centrality of (moderately) weighted networks and third a logarithmic variant of the eigenvector quantization in weighted networks.

Katz centrality

Apart from the eigenvector centrality (3), one often considers the Katz centrality which is defined as

$${x}_{i}^{(K)}=\alpha \sum _{j}\,{M}_{ij}{x}_{j}^{(K)}+\beta ,$$ (6)

where 0 < α < (λ 1 )−1 and β > 0. The main difference between the Katz and the eigenvector centralities is the second constant term which gives any vertex (even those without any in-link) a score β in addition to the contribution from the centrality of its neighbours. We would like to show the possibility to observe an approximate eigenvector-quantization-like picture for the Katz centrality. In fact, the histograms of single-cycle networks Fig. 4(b) still show the level structure, but as long as α decreases, these levels are shifted from integer values and split up into sublevels.

Figure 4 Left: histograms of the Katz centrality for single-cycle networks. The histograms show the number of nodes with a given value of \({x}^{(K)}/{x}_{{\rm{\min }}}^{(K)}\) for different values of α ∈ {1, 0.95, 0.9, 0.85} and β = 1. For α → 1, the same quantization can be seen. For α ≠ 1, while the level structure is still observed, these levels are no longer at integer values and some of them can be split up into sub-levels. Right: histograms of weighted networks. The single-cycle network from Fig. 1 right has been subjected to a weighting of its links. The weights are drawn from the interval [1 − ε, 1] uniformly at random. The histograms show the number of nodes with a given value of x/x min for various values of ε ∈ {0, 0.1, 0.3, 0.5}. For ε = 0 the pure quantization can be seen. For ε > 0 the signature deteriorates. Full size image

Slightly weighted networks

We consider a moderately weighted network, by drawing the weights of the links from the interval [1 − ε, 1] uniformly at random, where ε is a noise parameter. When ε = 0 we recover the unweighted case. We investigate to which extend it is possible to distinguish the single-cycle networks from the multicycle networks. The noise parameter ε is varied from zero, the baseline model, to one, the fully random model.

The histograms of single-cycle networks still show the quantization peaks as one departs from the baseline scenario, but as ε increases, these peaks broaden and gradually lose their equidistance. For the single-cycle network in Fig. 4 left, we have looked at the histograms of x of weighted networks for various values of ε. The examples are seen in Fig. 4 right.

Example of the weighted network — Leslie model of population aging

The Leslie model is used in population ecology to understand age distributions by means of a linear reproduction dynamics. Even though it is mathematically trivial, it is relevant ecologically45. For reproductions in ecosystems, it is important that cycles are present. The life-cycles of animals that undergo metamorphosis and significant transformations are examples of these reproductive cycles. Salmons, for instance, pass through various stages that involve the adaptation from sweet water to saltwater environments, and large migrations46. If at any of these stages there is a bottleneck, e.g. closed migration routes, extensive fishery, or a major ecological catastrophe, then the reproduction cycle is broken and the entire population is at risk of extinction.

We now describe how the theorem extends to the case of weighted links by means of the Leslie model. The population is divided into N stages, with population sizes x 1 (t), …, x N (t) at time t. The population at time t + 1 is obtained from the population at t by (see also Fig. 5)

$$\begin{array}{rcl}{x}_{1}(t+1) & = & \sum _{j > 1}\,{f}_{j}{x}_{j}(t)\\ {x}_{j}(t+1) & = & {p}_{j-1}{x}_{j-1}(t)\,{\rm{for}}\,j > 1.\end{array}$$ (7)

Figure 5 Illustration of the transitions in the Leslie model. There are N stages in the population (circles), each characterized by a population size x i . A percentage p i ≤ 1 of the population in the ith stage survives to the next stage p i+1 , while the number of new individuals born into stage 1 for a given individual at stage i is f i . Full size image

The first equation represents reproduction, where f j is the number of offspring per individual at stage j. The second set of equations determines the survival (or aging) process, where p j is the fraction of individuals that survive to the next stage. In matrix form, x(t + 1) = Lx(t), the components are given by L ij = f j δ i,1 + p j−1 δ i−1,j , where δ ij is the Kronecker delta. Here L is called the Leslie matrix and (7) is a discrete version of (1).

Since the Leslie matrix L is weighted, one cannot use the quantization criterion directly. However, often one is interested to know whether there are one or more cycles: certain species reproduce only in a certain stage of their life. This means that there is only one non-zero f i , and f j = 0, for j ≠ i. Let us also consider the case where all survival percentages are similar, such that p i ~ p, ∀i.

Without loss of generality, we choose the reproductive stage to be the last one, i = N. If this is not the case, we focus only on the part of the system, which creates the reproduction cycle; the remaining stages do not feed back into the cycle and thus don’t contribute to the reproduction. For long times, the population is determined by the eigenvector equation, Lx = λ 1 x. In this case, it is possible to calculate λ 1 and the corresponding eigenvector, \({\lambda }_{1}=p{(f/p)}^{\frac{1}{N}}\), and

$$x=({(\frac{f}{p})}^{\frac{N-1}{N}},{(\frac{f}{p})}^{\frac{N-2}{N}},\ldots ,{(\frac{f}{p})}^{\frac{1}{N}},1).$$

The logarithm of the population vector then reads

$$\log \,x=\frac{1}{N}\,\log \,\frac{f}{p}\cdot (N-1,N-2,\ldots ,1,0),$$ (8)