Applicability of spatial early warning signals to complex network dynamics

TL;DR

This paper explores how spatial early warning signals can predict tipping points in complex networks, finding that some signals work better than others.

Contribution

The study evaluates six spatial early warning signals on various networks, revealing their performance differences and reliability.

Findings

Coefficient of variation and spatial skewness outperform other early warning signals in many scenarios.

Spatial EWSs are more reliable on complex networks than on square lattices.

Performance of EWSs varies depending on the tipping scenario.

Abstract

Early warning signals (EWSs) for complex dynamical systems aim to anticipate tipping points before they occur. While signals computed from time-series data, such as temporal variance, are useful for this task, they are costly to obtain in practice because they need many samples over time to calculate. Spatial EWSs use just a single sample per spatial location and aggregate the samples over space rather than time to try to mitigate this limitation. However, although many complex systems in nature and society form diverse networks, the performance of spatial EWSs is mostly unknown for general networks because the vast majority of studies of spatial EWSs have been on regular lattice networks. Therefore, we have carried out a comprehensive investigation of six major spatial EWSs on various networks. We find that the winning EWS depends on tipping scenarios, although the coefficient of variation and spatial skewness tend to outperform alternative EWSs. We also find that spatial EWSs behave in a drastically different manner between the square lattice and complex networks and tend to be more reliable for the latter than the former. The present results encourage further studies of spatial EWSs on complex networks.

Figures

Overview of simulations and classification of EWSs. (a) Ascending sequence of simulations of the coupled double-well dynamics on a Barabási-Albert network with 100 nodes. Each grey line shows the xi* values at a node. When u is small, the xi* are all near their initial value, shown by the green line. As u increases, each xi* tends to increase but stays small until the first tipping point, annotated in red. After the first tipping point, at least some nodes move to a qualitatively different state. The home range of the control parameter is defined to be from the first control parameter value, here u=0, to the control parameter value immediately before the first tipping point. The simulation range contains the home range and encompasses all sampled values of the control parameter. (b) Descending sequence of simulations of the same dynamics on the same network. (c) Three hypothetical EWSs (brown lines) for the ascending simulations shown in (a). We classify an EWS by comparing linear regressions of samples of the EWS taken far from (line a, blue) and near (line b, red) the first tipping point in the home range. We classify the EWS as ‘accelerating’, ‘reversing’ or ‘unsuccessful‘ according to the criteria shown. (d) Three hypothetical EWS for the descending simulations are shown in (b).

Node states and EWSs as a function of D for the coupled double-well dynamics on the drug interaction network. (a) Ascending simulations. (b) Descending simulations.

Sign-corrected Kendall's τ values, τ′, for four dynamics models, two control parameters (i.e. D, shown by the circles, and u, shown by the crosses) and 35 networks. For each set of simulations, IM is shown in orange, CV in blue, g1′ in green and g2 in magenta. (a) Ascending simulations. (b) Descending simulations.

Classification of EWSs for four dynamic models on 35 networks. Solid bars indicate the fraction of networks for which the EWS is successful (either accelerating or reversing category). (a) Ascending simulations. (b) Descending simulations

Spatial coefficient of variation (CV) and spatial skewness ( g1′) as a control parameter varies towards the tipping point. Each grey line shows xi*. The spread of the xi* obtained at four control parameter values is shown in black and labelled with the values of CV (blue) and g1′ (green) at the control parameter values. Classification results are provided in the corresponding coloured text. The direction of the simulation sequence is indicated by the arrow beneath the horizontal axis. We used the drug interaction network. (a) and (b): coupled double-well, with u as the control parameter. (c) and (d): coupled double-well, with D as the control parameter. (e) and (f): SIS, with D as the control parameter.

IM and the two quantities defining IM for the simulations shown in figure 5. In each panel, we show the value of IM (orange line), its numerator (blue line) and its denominator (red line). The direction of the simulation sequence is indicated by the arrow beneath the horizontal axis. The classification result for IM is shown in orange text. (a) and (b): coupled double-well, with u as the control parameter. (c) and (d): coupled double-well, with D as the control parameter. (e) and (f): SIS, with D as the control parameter.

Equations

Full text

Introduction

Complex systems display tipping points when there exists some environmental threshold beyond which the system enters a qualitatively different regime [1]. For example, tropical woodland ecosystems may collapse to a relatively barren state as rainfall decreases across a critical threshold [2]. As another example, a novel communicable disease may start to rapidly spread in a population when some environmental conditions are met [3]. Tipping points are in general difficult to anticipate because small changes in driver variables can have markedly different effects on the state of the system [1]. However, a variety of systems display characteristic behaviours in the proximity of a tipping point, and such behaviours have been exploited for developing several early warning signals (EWS) that can anticipate the onset of a tipping point [4–7].

Dynamical systems near a tipping point recover more slowly from a disturbance than those far from a tipping point. This phenomenon, called critical slowing down, leads to increased autocorrelation and variance in time-series data, which are typical EWSs [4]. In fact, to calculate these EWSs, many samples from the same element in the system are required in each environmental condition (e.g. a control parameter value in the case of modelling studies) and over several environmental conditions (i.e. some far from a tipping point and others nearer to it) [6]. For example, when samples independently obey an identical normal distribution, emulating one environmental condition, the sample standard deviation, which is a typical EWS, has a standard deviation proportional to $[eqn]$ , where $[eqn]$ is the number of samples [8–11]. Therefore, ideally, one wants to secure $[eqn]$ samples or more to reliably estimate the sample standard deviation. However, in practice, it is often too costly to collect so many samples per environmental condition [12,13], potentially contributing to the lack of consistent EWSs in empirical systems [14,15]. Furthermore, if $[eqn]$ is large, the environment may drift to a different state in the middle of collecting $[eqn]$ samples in the field or experiment. If this is the case, the use of the EWS computed from the $[eqn]$ samples is compromised because the EWS reflects a range of environmental conditions rather than a single one.

Spatial EWSs seek to mitigate this limitation of ‘temporal’ EWSs by measuring, in each environmental condition, a single sample from many different elements constituting a complex system, rather than obtaining many samples from one element (or multiple elements) in the system [7,13]. We define spatial EWSs as requiring just one sample per element, i.e. $[eqn]$ . Several proposed spatial EWSs are spatial analogues of temporal EWSs, such as spatial correlation [13], spatial variance and skewness [16], the power spectrum of a state variable of a spatially extended system [17] and recovery length (i.e. as opposed to recovery time) [18]. Other spatial EWSs have no temporal analogue, such as the distribution of patch sizes in patchy environments [19].

EWSs have been probably most vigorously studied for ecological dynamics, many of which take place in physical space. Presumably for this reason, most studies of spatial EWSs have been carried out on spatial regular, grid-like networks modelling two-dimensional ecological landscapes, such as the square lattice with or without periodic boundary conditions [7,19–24] and partial differential equations involving space and time [16,17,20,25–30]. Such a network was also used in a study of EWSs for deforestation transitions [2]. However, many other empirical complex systems for which prediction of tipping points is desired have more complex network structure than regular lattices or two-dimensional continuous planes. Examples include epidemic spreading in human and animal populations [3], progression of diseases in general [5] and symptoms of mental disorders in particular [31,32], and inter-specific population dynamics among animals and plants [33]. Furthermore, even if ecological dynamics occurs in a two-dimensional terrain, habitats may be irregularly distributed and heterogeneous [34,35] such that the underlying network is spatial but heterogeneous. Therefore, empirical studies of spatial EWSs in ecological systems [18,19,27,36–44] may be better justified if spatial EWSs are shown to be valid on heterogeneous networks rather than on regular lattices. However, spatial EWSs have been rarely studied beyond on regular lattices, while notable exceptions exist for complex dynamics models coupling epidemic and opinion dynamics [45,46].

In summary, despite the need, whether or not and which spatial EWSs perform well on heterogeneous networks is largely unknown. In the present study, we comprehensively investigate the performance of six major spatial EWSs on complex networks, compared across dynamics models, environmental parameters, how the tipping point is approached and networks. We also provide mechanistic understanding of why they work or do not work depending on the situation and recommended practices based on our numerical results.

Methods

Dynamics

2.1.

We used the following four stochastic dynamics models on networks: a coupled double-well model, a model of mutualistic species interactions, a susceptible–infectious–susceptible (SIS) model and a gene-regulatory model. For the coupled double-well, mutualistic species and gene-regulatory dynamics models, we consider two control parameters: the strength of coupling between nodes, $[eqn]$ and a stress parameter, $[eqn]$ , which can exert negative ( $[eqn]$ ) or positive ( $[eqn]$ ) stress uniformly on all nodes. For the SIS model, we only consider $[eqn]$ , which is more conventionally known as the infection rate, as the sole control parameter because the concept of a uniform stress is not realistic for the SIS model. The matrix $[eqn]$ is the adjacency matrix of the network. We assume that the network is connected, undirected and unweighted. Each dynamics is simulated with Gaussian noise $[eqn]$ with noise strength $[eqn]$ . The noise applied to different nodes is assumed to be independent of each other.

A coupled double-well model on networks is given by [47]:

$$[eqn]$$

where $[eqn]$ is the dynamical state of the $[eqn]$ th node and represents a numeric attribute, such as the species population size or amount of tree cover; $[eqn]$ are parameters which, in the absence of noise and coupling, set the positions of the equilibria and $[eqn]$ is the number of nodes. The coupled double-well dynamics has been used for modelling various phenomena, including human social movements [48], interacting biological species [49] and connected climate regions [2]. In the absence of coupling and noise, there are two stable equilibria: a lower equilibrium at $[eqn]$ and an upper equilibrium at $[eqn]$ . We set $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ (if the control parameter is $[eqn]$ ), $[eqn]$ (if the control parameter is $[eqn]$ ) and $[eqn]$ . Exceptions to parameter values are provided in the electronic supplementary material. We initialize this model in either the lower state with $[eqn]$ $[eqn]$ or the upper state with $[eqn]$ $[eqn]$ . We consider an $[eqn]$ th node to be in the lower state if $[eqn]$ and in the upper state otherwise.

A model of mutualistic species dynamics is given by [50]:

$$[eqn]$$

where $[eqn]$ represents the abundance of the $[eqn]$ th species, $[eqn]$ is a constant incoming migration rate, $[eqn]$ is the carrying capacity, $[eqn]$ is the Allee constant and $[eqn]$ , $[eqn]$ and $[eqn]$ moderate the effect of the interaction term $[eqn]$ . By following [50], we set $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ , $[eqn]$ and $[eqn]$ . We use $[eqn]$ (if the control parameter is $[eqn]$ ), $[eqn]$ (if the control parameter is $[eqn]$ ), and $[eqn]$ . In the absence of coupling and dynamical noise, this model has a stable lower state with $[eqn]$ and a stable upper state with $[eqn]$ . We initialize this model in the upper state with $[eqn]$ $[eqn]$ , which is an arbitrary large value representing an extant population [50]. We start the dynamics from the upper state because a common interest in ecology is loss of resilience in the current species assemblage, which is modelled by collapse from the upper state [4]. We consider that any $[eqn]$ and any $[eqn]$ are in the upper and lower state, respectively. To prevent $[eqn]$ for any node $[eqn]$ and time $[eqn]$ , which is not physical for this model, we set $[eqn]$ whenever our quadrature algorithm produced $[eqn]$ during simulations. We also use the same procedure to prevent $[eqn]$ for the following two models.

An SIS model on networks in the stochastic differential equation form is given by [3]:

$$[eqn]$$

The node state $[eqn]$ represents the probability that the $[eqn]$ th node is infectious (the $[eqn]$ th node is susceptible with probability $[eqn]$ ); $[eqn]$ is the infection rate and $[eqn]$ is the recovery rate. We use $[eqn]$ and $[eqn]$ . In the absence of noise, there is a disease-free equilibrium with $[eqn]$ $[eqn]$ , which always exists and is stable when $[eqn]$ is below an epidemic threshold value, and an endemic equilibrium in which $[eqn]$ $[eqn]$ , which exists and is stable when $[eqn]$ is large enough [3]. In the presence of noise, some $[eqn]$ are expected at any value of $[eqn]$ . We simulate this model beginning in either the lower (i.e. almost disease-free) state with $[eqn]$ $[eqn]$ or the upper (i.e. endemic) state with $[eqn]$ . We consider an $[eqn]$ th node to be in the lower state if $[eqn]$ and in the upper state otherwise.

A model of gene-regulatory dynamics on networks is given by [50]:

$$[eqn]$$

where $[eqn]$ represents the expression level of the $[eqn]$ th gene, $[eqn]$ and $[eqn]$ characterize the behaviour of the $[eqn]$ th gene in isolation and $[eqn]$ controls the interaction of the $[eqn]$ th and $[eqn]$ th genes. Following [50], we set $[eqn]$ , $[eqn]$ and $[eqn]$ . We use $[eqn]$ (if the control parameter is $[eqn]$ ), $[eqn]$ (if the control parameter is $[eqn]$ ) and $[eqn]$ . In the absence of noise, this model has an equilibrium at $[eqn]$ $[eqn]$ , which is stable when $[eqn]$ or $[eqn]$ is small enough and represents the inactive state, and an active state with $[eqn]$ $[eqn]$ , which exists when $[eqn]$ or $[eqn]$ is sufficiently large. We simulate this model from the upper state with $[eqn]$ $[eqn]$ because one is often interested in modelling the loss of resilience of the active state [50]. We use the same criteria to define the lower and upper states for the gene-regulatory dynamics as we did for the SIS dynamics.

Networks

2.2.

We used 30 empirical networks and five synthetic networks built with different generative models. These networks vary in terms of the number of nodes and edges, the heterogeneity of the degree distribution and community structure. Separately, we also analysed a square lattice for comparison purposes. We coerce each network to be undirected and unweighted if it is not and use only the largest connected component. Details of the individual networks are found in the electronic supplementary material, section S1.

Early warning signals

2.3.

We compare six types of EWS that require only $[eqn]$ sample of $[eqn]$ for each node $[eqn]$ at any given control parameter value, i.e. Moran’s $[eqn]$ , the spatial standard deviation, spatial variance, spatial coefficient of variation (CV), spatial skewness and spatial kurtosis. We compute the EWSs using the equilibrium $[eqn]$ values defined below. Furthermore, the information required by these EWSs is modest: all but Moran’s $[eqn]$ need only the $[eqn]$ values; Moran’s $[eqn]$ needs the $[eqn]$ values and the adjacency matrix of the network. We chose these six EWSs because spatial variance is a straightforward variant of the spatial standard deviation and because, apart from the EWSs specific to ecological populations, the other five EWSs are commonly used in the literature.

Moran’s $[eqn]$ is defined as follows [51,52]:

$$[eqn]$$

where $[eqn]$ and $[eqn]$ . Moran’s $[eqn]$ is a measure of spatial correlation because it quantifies the extent to which neighbouring sampling sites on the same surface or object (i.e. nodes in the case of networks [52]) have similar states [51]. Specifically, $[eqn]$ is the ratio of the normalized cross-products, or covariance, of the node states, $[eqn]$ , to their total variance, $[eqn]$ . Moran’s $[eqn]$ is similar to Pearson’s correlation coefficient in that it is close to 1 when neighbouring $[eqn]$ have similar values, close to $[eqn]$ when they are dissimilar and close to 0 when they are not correlated. However, $[eqn]$ can be less than $[eqn]$ or more than $[eqn]$ [51].

We refer to the sample standard deviation of the node states as the spatial standard deviation, denoted by $[eqn]$ , i.e.:

$$[eqn]$$

Quantity $[eqn]$ is the unbiased estimate of the population standard deviation and is often used as an EWS in spatially extended systems [26,39,40]. We also use the sample variance, $[eqn]$ , where

$$[eqn]$$

is the $[eqn]$ th central moment of $[eqn]$ , and the CV. The CV is the normalized variant of $[eqn]$ defined by

$$[eqn]$$

We note that the CV is often used as a spatial EWS [18,36,40,41].

Skewness and kurtosis have been proposed as EWS for time-series data [6], and both have also been used as spatial EWS [16,26,29,39]. Skewness and kurtosis are the scaled third and fourth central moments, respectively, of the probability distribution of a random variable $[eqn]$ . Sample skewness is defined as

$$[eqn]$$

and quantifies the extent to which extreme values tend to appear to the right (positive) or left (negative) side of the mean and approaches 0 as $[eqn]$ for a symmetric distribution. Skewness may increase (i.e. become more positive) or decrease (i.e. become more negative) as a system approaches a tipping point [7]. Therefore, to ensure that the desired behaviour of $[eqn]$ as a system approaches a tipping point is encoded into a more positive value, we use $[eqn]$ as the EWS, where $[eqn]$ for ascending simulations and $[eqn]$ for descending simulations (see §2.4 for the simulation protocol and §2.5 for a similar procedure involving Kendall’s $[eqn]$ ).

Sample kurtosis is defined as

$$[eqn]$$

and quantifies the magnitude of the extreme values (i.e. how large and how far from the mean) and approaches 3 as $[eqn]$ for a normal distribution. An adjustment $[eqn]$ can be used to define excess kurtosis with respect to a standard normal distribution, but we do not use that convention here. Because extreme values should become more common near a tipping point, a larger (i.e. more positive) kurtosis may indicate an approaching tipping point [23].

Simulations

2.4.

We performed simulations for each combination of dynamics model (i.e. double-well, mutualistic species, SIS or gene-regulatory), control parameter (i.e. $[eqn]$ or $[eqn]$ ), direction (ascending or descending; see below) and network (i.e. one of the 35 networks) and measured the performance of the six EWSs. We refer to a combination of a dynamics model, control parameter and direction as the simulation condition. As an example, we consider the coupled double-well dynamics with $[eqn]$ as the control parameter initialized in the lower state (corresponding to the ‘ascending’ simulations as explained below), which altogether specifies a simulation condition, on a Barabási-Albert network (figure 1a). We conduct the first simulation with $[eqn]$ , setting $[eqn]$ $[eqn]$ . We integrate equation (2.1) using the Euler–Maruyama method with $[eqn]$ for 50 time units. We consider $[eqn]$ as an equilibrated value under the dynamical noise, collect one sample of $[eqn]$ from each $[eqn]$ th node, and calculate the EWSs. Then, we increase $[eqn]$ by a small amount and perform the simulation again using the same initial state. The procedure when $[eqn]$ is the control parameter is similar.

Overview of simulations and classification of EWSs. (a) Ascending sequence of simulations of the coupled double-well dynamics on a Barabási-Albert network with 100 nodes. Each grey line shows the xi values at a node. When u is small, the xi* are all near their initial value, shown by the green line. As u increases, each xi* tends to increase but stays small until the first tipping point, annotated in red. After the first tipping point, at least some nodes move to a qualitatively different state. The home range of the control parameter is defined to be from the first control parameter value, here u=0, to the control parameter value immediately before the first tipping point. The simulation range contains the home range and encompasses all sampled values of the control parameter. (b) Descending sequence of simulations of the same dynamics on the same network. (c) Three hypothetical EWSs (brown lines) for the ascending simulations shown in (a). We classify an EWS by comparing linear regressions of samples of the EWS taken far from (line a, blue) and near (line b, red) the first tipping point in the home range. We classify the EWS as ‘accelerating’, ‘reversing’ or ‘unsuccessful‘ according to the criteria shown. (d) Three hypothetical EWS for the descending simulations are shown in (b).*

We iterate these steps to simulate the dynamics and calculate EWSs at $[eqn]$ evenly spaced control parameter values. In the present case, we use $[eqn]$ values of $[eqn]$ in the range $[eqn]$ , which we call the simulation range (figure 1a). This range of $[eqn]$ results in $[eqn]$ $[eqn]$ when $[eqn]$ is smaller than $[eqn]$ and at least some $[eqn]$ in a qualitatively different state, i.e. the upper state (i.e. $[eqn]$ ), when $[eqn]$ is larger.

At large $[eqn]$ , the statistics of $[eqn]$ of the nodes in the upper state are not useful for predicting which nodes will next make the transition from the lower to the upper state [53]. Therefore, we analyse each EWS only for the control parameter values for which all $[eqn]$ values are near its initial state (i.e. the lower state in the present case). We refer to this restricted parameter range as the home range. As shown in figure 1a, the home range is a subset of the simulation range. Note that the home range is specific to each combination of the simulation condition and network.

We ensured the following two properties in our simulations. First, the simulation range contains at least one tipping point such that the home range is well defined (i.e. the last control parameter value in the home range is just before the first tipping point). Second, the simulation range contains sufficiently many control parameter values near and far from the tipping point such that we can assess the performance of the EWSs as described in the following text. To ensure these two properties, we determined the simulation range by trial and error separately for each simulation condition and network. We show the simulation range for each simulation condition and network in the electronic supplementary material, section S2.

The sequence of simulations shown in figure 1a begins with each $[eqn]$ in the lower state and ends when at least some $[eqn]$ enter the upper state. We call such a sequence of simulations an ascending simulation. By contrast, in figure 1b, we initialize each $[eqn]$ in the upper state, i.e. $[eqn]$ $[eqn]$ , and keep reducing $[eqn]$ by a small amount for each of the $[eqn]$ simulations. We call this type of sequence of simulations a descending simulation.

Kendall’s τ

2.5.

As in prior work, we used the Kendall’s rank correlation, denoted by $[eqn]$ , as a performance measure of EWSs [7,13]. We computed $[eqn]$ between the control parameter and the EWS over a range of control parameter values near the tipping event, specifically, the latter half of the home range. We also computed a sign-adjusted Kendall’s $[eqn]$ value as follows. For dynamics simulated in an ascending sequence, a positive rank correlation indicates that the EWS grows large as the control parameter approaches a bifurcation. In this case, we set $[eqn]$ . For dynamics simulated in a descending sequence, a good EWS should grow large as the control parameter becomes smaller (more negative). In this case, we set $[eqn]$ . A larger $[eqn]$ value is better. Both $[eqn]$ and $[eqn]$ range between $[eqn]$ and $[eqn]$ .

Classification of early warning signals

2.6.

Kendall’s $[eqn]$ is a dominant performance measure for EWSs but with criticisms [54,55]. Our numerical simulations produced diverse behaviour of the different EWSs as we vary the control parameter, including the case in which the EWS decreases as we approach the tipping point. Given this situation, solely relying on Kendall’s $[eqn]$ would not generate useful comparison between EWSs. Therefore, we developed a classification scheme of EWS as follows.

Consider the example simulation sequence shown in figure 1a, in which we start with the lower state and gradually increase the control parameter, $[eqn]$ . Suppose that an EWS responds to the gradual increase in $[eqn]$ within the home range of $[eqn]$ as shown by the uppermost brown line in figure 1c. This EWS is noisy but remains low when $[eqn]$ . It then increases progressively rapidly as $[eqn]$ approaches the tipping point. This is an ideal case because small values of the EWS reliably indicate that the system is far from transition, whereas large values indicate that the first transition is nearby in terms of $[eqn]$ . We quantify the extent to which an EWS follows this pattern using the trend in the EWS value when it is far from the first transition (slope $[eqn]$ , blue line in the upper panel of figure 1c) and when it is near (slope $[eqn]$ , red line). In the upper panel in figure 1c, both $[eqn]$ and $[eqn]$ are positive and $[eqn]$ is substantially larger than $[eqn]$ . We classify an EWS that behaves in this fashion as successful and say that the EWS is accelerating. We will give the precise definition below.

Alternatively, an EWS may first tend to decrease as $[eqn]$ increases when $[eqn]$ is far from its tipping point, as in the middle panel of figure 1c. However, if the EWS steadily increases at larger $[eqn]$ values near the tipping point, this EWS behaviour also reliably indicates that the system is approaching an impending tipping point. Therefore, we also classify this behaviour as successful and say that the EWS is reversing.

Many other trends in EWS behaviour as a function of the control parameter are possible. For example, the EWS value may initially rise rapidly as $[eqn]$ increases far from a tipping point. Then, the EWS may level off or even decrease as $[eqn]$ further increases, approaching the tipping point, as in the lower panel of figure 1c. Such an EWS gives a false positive while the system is still far from the bifurcation and a false negative when it is close to the bifurcation. We classify such a behaviour, and any other pattern not covered by the accelerating and reversing categories, as unsuccessful.

We classify EWSs in the same manner when we start from an upper state and gradually decrease the control parameter, as we show in figure 1d. The lower panel of figure 1d shows a different case of failure (i.e. neither category) from that shown in the lower panel of figure 1c just for demonstration.

To compute $[eqn]$ , we run an ordinary linear least square regression on the first five values of the control parameter from the home range, in which the independent variable is the control parameter and the dependent variable is the EWS. We compute $[eqn]$ in the same fashion but using the last five values of the control parameter in the home range. If $[eqn]$ and $[eqn]$ are positive and $[eqn]$ , where $[eqn]$ , we say that the EWS is accelerating. We set $[eqn]$ . If $[eqn]$ is negative, $[eqn]$ is positive and $[eqn]$ , then we say that the EWS is reversing. We use $[eqn]$ to capture trends in an EWS that have increased markedly with respect to an initial negative trend without being too strict.

Results

We ran numerical simulations on four dynamical system models on networks, i.e. coupled double-well, mutualistic species, SIS and gene-regulatory models. We performed sequences of simulations of each model across a range of parameter values, forcing a bifurcation, on 35 empirical and model networks. We then computed the six spatial early warning signals for each simulation sequence. To assess the quality of each EWS, we computed $[eqn]$ , the sign-corrected version of Kendall’s rank correlation. We further classified each EWS as accelerating, reversing or unsuccessful based on the extent to which it showed desirable warning signal behaviour far from and near to the tipping point. We found that the overall trends of the performance results were similar among the three EWSs based on the spatial variance, i.e. $[eqn]$ , $[eqn]$ and $[eqn]$ , and that $[eqn]$ outperformed $[eqn]$ and $[eqn]$ in most simulations (see the electronic supplementary material, S3, S4 for the results). Therefore, we do not show the results for $[eqn]$ or $[eqn]$ in the following text unless we state otherwise.

Examples

3.1.

As an example, let us consider the coupled double-well dynamics (see equation (2.1)) on the drug interaction network. We use the coupling strength, $[eqn]$ , as the control parameter, which we gradually increase starting from zero (i.e. ascending simulations). The grey lines in figure 2a represent the $[eqn]$ values for all the nodes as a function of $[eqn]$ . When $[eqn]$ is small, all the nodes are in their lower state (i.e. $[eqn]$ ). As $[eqn]$ gradually increases but is still smaller than the tipping point, each $[eqn]$ becomes larger but still remains in the lower state. The first transition of any node to the upper state occurs around $[eqn]$ , and progressively larger values of $[eqn]$ result in the transition of more nodes.

Node states and EWSs as a function of D for the coupled double-well dynamics on the drug interaction network. (a) Ascending simulations. (b) Descending simulations.

The orange line in figure 2a indicates that $[eqn]$ initially increases rapidly, then levels off and even decreases as $[eqn]$ approaches the tipping point. This pattern of $[eqn]$ values is not desirable as an EWS and is reflected in the negative sign-adjusted Kendall’s rank correlation value ( $[eqn]$ ); our algorithm classifies $[eqn]$ as ‘unsuccessful.’

By contrast, $[eqn]$ , shown by the blue line, grows in an accelerating fashion as $[eqn]$ approaches the tipping point, yielding $[eqn]$ and classification as ‘accelerating’, one of our two successful categories. Both $[eqn]$ (green) and $[eqn]$ (magenta) behave similarly to, while more noisily than, $[eqn]$ , except that $[eqn]$ and $[eqn]$ rapidly increase initially as $[eqn]$ increases. Although both $[eqn]$ and $[eqn]$ are desirable EWSs in terms of $[eqn]$ , with $[eqn]$ , our classification scheme classifies them as ‘unsuccessful’ owing to their initial rapid increase. However, if the simulation starts from a larger value of $[eqn]$ while the last value of $[eqn]$ remains the same, we can remove the part of the curve in figure 2a where $[eqn]$ or $[eqn]$ rapidly increases. Then, $[eqn]$ and $[eqn]$ would be classified to be successful. See the electronic supplementary material, S4, for results supporting this phenomenon across different networks. Given such a scenario, we advise to refer to both $[eqn]$ and our classification scheme in general.

The performance of each EWS depends on the simulation condition (i.e. combination of a dynamics, control parameter and ascending versus descending simulation direction). In figure 2b, we show each $[eqn]$ when we initialize $[eqn]$ in the upper state and gradually decrease $[eqn]$ from $[eqn]$ (i.e. descending simulations); the dynamics model and network are the same as that used in figure 2a. In this case, $[eqn]$ initially decreases when $[eqn]$ decreases and is still far from the tipping point. Then, $[eqn]$ markedly increases as $[eqn]$ further decreases and approaches the tipping point. Overall, the trend of $[eqn]$ values is desirable; we obtain $[eqn]$ and classify $[eqn]$ as ‘reversing’, one of the two successful categories.

The other EWSs also perform differently from how they did in figure 2a. In figure 2b, $[eqn]$ performs poorly in terms of $[eqn]$ ; its overall trend is largely linear and monotonically decreasing (rather than increasing) towards the tipping point, yielding $[eqn]$ . We classify $[eqn]$ as ‘reversing’ because $[eqn]$ increases near the bifurcation. Finally, $[eqn]$ behaves nearly ideally ( $[eqn]$ , classified as ‘accelerating’), whereas $[eqn]$ behaves almost conversely ( $[eqn]$ , classified as ‘unsuccessful’).

In summary, the performance of the different EWSs can substantially vary depending on the direction of the simulations (i.e. ascending versus descending simulations) although the dynamics model and the network are the same. Therefore, we expect that there are various situations in which one EWS may work better than another and vice versa, which we investigate in the following sections.

Variability in early warning signals performance over simulation conditions

3.2.

To assess the variation in performance of the different EWSs across conditions, we carried out simulations with each simulation condition on 35 networks.

As an initial analysis, we show the $[eqn]$ values for each EWS and simulation condition in figure 3. Note that we omitted six simulation conditions as being unrealistic (see §2 ). Therefore, the results for these simulation conditions are missing in figure 3 (i.e. mutualistic species and gene-regulatory in figure 3a and variation in $[eqn]$ for the SIS). We find that $[eqn]$ sometimes performs well (i.e. shown by orange markers near $[eqn]$ ) and sometimes poorly ( $[eqn]$ much smaller than $[eqn]$ , including near $[eqn]$ ). There are many intermediate values of $[eqn]$ , particularly for ascending simulations (figure 3a). On the other hand, $[eqn]$ performs remarkably well when the control parameter is $[eqn]$ (i.e. stress) and under some other simulation conditions in which the control parameter is $[eqn]$ . Skewness $[eqn]$ generally performs well on the coupled double-well and mutualistic species dynamics, and its performance is mixed on the SIS and gene-regulatory dynamics; $[eqn]$ performs best on ascending simulations of the coupled double-well dynamics and has mixed performance otherwise.

Sign-corrected Kendall's τ values, τ′, for four dynamics models, two control parameters (i.e. D, shown by the circles, and u, shown by the crosses) and 35 networks. For each set of simulations, IM is shown in orange, CV in blue, g1′ in green and g2 in magenta. (a) Ascending simulations. (b) Descending simulations.

We applied our classification procedure to further quantify the performance of the different EWSs. The solid colour bars in figure 4 represent the proportion of networks out of the 35 networks with successful EWSs (i.e. classified as either accelerating or reversing) under each simulation condition. The figure supports the observation made in figure 3 that the performance of an EWS depends on the simulation condition. For example, $[eqn]$ is successful for the largest proportion of networks on average among all the EWSs, including $[eqn]$ and $[eqn]$ (see the electronic supplementary material, S4 for the results). Instead, $[eqn]$ is the best for the coupled double-well and mutualistic species dynamics in descending simulations regardless of the control parameter. Although figure 4 indicates that $[eqn]$ and $[eqn]$ perform poorly in ascending simulations of the coupled double-well dynamics when the control parameter is $[eqn]$ , the performance considerably improves if we start simulations from a larger $[eqn]$ value (see the electronic supplementary material, S4 for the results). Moran’s $[eqn]$ , on the other hand, is reliable for the SIS and gene-regulatory dynamics in the descending simulations.

Classification of EWSs for four dynamic models on 35 networks. Solid bars indicate the fraction of networks for which the EWS is successful (either accelerating or reversing category). (a) Ascending simulations. (b) Descending simulations

Figure 4 suggests that $[eqn]$ and $[eqn]$ perform better than $[eqn]$ and $[eqn]$ overall. To quantitatively investigate the relative performances of the EWSs including $[eqn]$ and $[eqn]$ , we assign each EWS a score based on its rank: the best EWS among the six EWSs for a particular simulation condition received 5 points, the second best a 4 and so forth and the worst a 0. In the case of a tie, the tied EWSs receive the average score (e.g. if two EWSs are tied for best, both receive 4.5 points). Summing this score across the 10 simulation conditions confirms that $[eqn]$ ( $[eqn]$ points) and $[eqn]$ ( $[eqn]$ points) are almost equally the best, followed by $[eqn]$ ( $[eqn]$ ), $[eqn]$ ( $[eqn]$ ), $[eqn]$ ( $[eqn]$ ) and $[eqn]$ ( $[eqn]$ ). Overall, $[eqn]$ was categorized as successful in most cases (76.3%), followed by $[eqn]$ (69.4%), $[eqn]$ (54.3%), $[eqn]$ (53.7%), then $[eqn]$ (48.0%) and $[eqn]$ (42.3%), as shown in the electronic supplementary material, S4. In terms of $[eqn]$ , it turned out that $[eqn]$ was the best, $[eqn]$ was the second best and the other four EWSs performed notably worse than $[eqn]$ and $[eqn]$ (see the electronic supplementary material, S3 for the results). Therefore, it appears that $[eqn]$ and $[eqn]$ are the most reliable spatial EWSs among the six EWSs.

In §1, we pointed out that most studies of spatial EWSs were carried out on the square lattice or its continuous variants. Therefore, we examined the performances of the spatial EWSs on the square lattice with $[eqn]$ nodes and periodic boundary conditions. We find that results for the square lattice are substantially different from those for the 35 networks. Specifically, our algorithm classifies $[eqn]$ as a successful EWS for the square lattice in 8 out of 10 simulation conditions, $[eqn]$ in 5 out of 10 conditions, $[eqn]$ in 2 out of 10 conditions and $[eqn]$ in 1 out of 10 conditions. The failure of $[eqn]$ is in stark contrast with the case of heterogeneous networks. We hypothesized that this last result is because all nodes are structurally equivalent to each other in the square lattice, such that $[eqn]$ values are statistically the same for all $[eqn]$ , yielding the lack of strong asymmetry in the distribution of $[eqn]$ . We have verified that, under the simulation conditions in which $[eqn]$ is unsuccessful, $[eqn]$ is centred around $[eqn]$ as the control parameter varies (see the electronic supplementary material, S5). Furthermore, except $[eqn]$ , spatial EWSs on the square lattice are typically noisy, leading to smaller $[eqn]$ values even when an EWS is classified as successful. For $[eqn]$ , the average $[eqn]$ value is $[eqn]$ when it is successful on the square lattice and $[eqn]$ when it is successful on all other networks. For $[eqn]$ , the average $[eqn]$ value is $[eqn]$ when it is successful on the square lattice and $[eqn]$ when it is successful on all other networks. We show examples of the noisiness of EWSs on the square lattice in the electronic supplementary material, S6. In summary, we conclude that what we know about spatial EWSs on the square lattice does not translate to the case of heterogeneous networks, which further strengthens the need for the present study.

Mechanisms of variable performance of early warning signals

3.3.

We have shown that different EWSs perform better than others in different simulation conditions. Focusing on heterogeneous networks, we explore mechanisms behind this observation in this section. Because $[eqn]$ was the worst performer in all our comparisons, we do not investigate it further.

Coefficient of variation

3.3.1.

To examine why $[eqn]$ performs well in some simulation conditions and not in others, let us consider again the coupled double-well dynamics on the drug interaction network used for the demonstration in figure 2. In figure 5a, we gradually increase $[eqn]$ and observe $[eqn]$ and $[eqn]$ (and $[eqn]$ ). In this case, $[eqn]$ increases in an accelerating fashion as $[eqn]$ approaches the tipping point, successfully anticipating the saddle-node bifurcation, which the $[eqn]$ values shown in blue in the figure at four values of $[eqn]$ indicate. Qualitatively the same behaviour occurs for all 35 networks. Figure 5a also indicates that the accelerated increase in $[eqn]$ is caused by the behaviour of the $[eqn]$ values of high-degree nodes (shown by the grey lines with larger $[eqn]$ ). These nodes receive more input from their neighbours and thus approach a bifurcation earlier than low-degree nodes do [53]. These $[eqn]$ values thus separate from those of the smaller degree nodes near the bifurcation, causing $[eqn]$ to increase.

Spatial coefficient of variation (CV) and spatial skewness ( g1′) as a control parameter varies towards the tipping point. Each grey line shows xi. The spread of the xi* obtained at four control parameter values is shown in black and labelled with the values of CV (blue) and g1′ (green) at the control parameter values. Classification results are provided in the corresponding coloured text. The direction of the simulation sequence is indicated by the arrow beneath the horizontal axis. We used the drug interaction network. (a) and (b): coupled double-well, with u as the control parameter. (c) and (d): coupled double-well, with D as the control parameter. (e) and (f): SIS, with D as the control parameter.*

In figure 5b, we show numerical results under the same simulation condition as in figure 5a except that the direction is reversed (i.e. descending simulations). In this case, small-degree nodes separate from the rest as the bifurcation is approached. Nevertheless, their effect on $[eqn]$ is the same as in the case of figure 5a. In other words, the value of $[eqn]$ grows owing to a progressive deviation of the smallest $[eqn]$ values at small-degree nodes as $[eqn]$ decreases towards the saddle-node bifurcation.

When the control parameter is $[eqn]$ , the behaviour of $[eqn]$ is not equivalent in both directions. In the ascending simulations, the $[eqn]$ values of large-degree nodes separate from the remainder (see figure 5c; this is the same simulation condition as in figure 2a), which is similar to when the control parameter is $[eqn]$ (see figure 5a). Therefore, $[eqn]$ increases in an accelerated manner, successfully anticipating the bifurcation. However, in descending simulations, the $[eqn]$ values of the different nodes come closer together as $[eqn]$ decreases towards the bifurcation (see figure 5d; also see figure 2b). This result is opposite to the results when the control parameter is $[eqn]$ (figure 5b). The reason for this behaviour is in that the input from the other nodes to the $[eqn]$ th node, given by $[eqn]$ , is roughly proportional to the degree of the $[eqn]$ th node (i.e. the number of $[eqn]$ values for which $[eqn]$ ). Owing to the multiplicative factor $[eqn]$ , as $[eqn]$ becomes smaller, the nodes receive smaller and more homogeneous amounts of input from their neighbours. Therefore, the $[eqn]$ values come closer, reducing $[eqn]$ . Although $[eqn]$ is still categorized as successful in figure 5d, it is because $[eqn]$ starts to increase very near the bifurcation (i.e. roughly between the two leftmost vertical thick lines in the figure) as $[eqn]$ decreases. The $[eqn]$ values behave similarly as $[eqn]$ gradually increases or decreases in the SIS dynamics, which show transcritical bifurcations (see figure 5e,f). However, whether or not $[eqn]$ is successful in figure 5e,f does not coincide with the results for the coupled double-well dynamics shown in figure 5c,d, respectively. This discrepancy is probably because $[eqn]$ is a quantity normalized by the mean of $[eqn]$ .

Spatial skewness

3.3.2.

The causes of the favourable behaviour of $[eqn]$ as an EWS are related to those of $[eqn]$ as follows.

First, $[eqn]$ and $[eqn]$ are successful in figure 5a,b for the same reason. EWS $[eqn]$ is also successful in figure 5e for the same reason although the increasing trend near the tipping point is weak. Specifically, as the tipping point is approached, $[eqn]$ of nodes closer to the tipping (i.e. high-degree nodes in figure 5a,e, and low-degree nodes in figure 5b) more rapidly deviate from the remainder, i.e. becoming larger in figure 5a,e, and smaller in figure 5b. Then, the overall variability of $[eqn]$ grows, captured as an increase in $[eqn]$ as well as $[eqn]$ . In this manner, skew increases in ascending simulations at least near the tipping point (i.e. figure 5a,e; increasing $[eqn]$ ), driven by the large-degree nodes, whereas skew decreases in descending simulations (i.e. figure 5b; this also increases $[eqn]$ ) as the low-degree nodes separate from the other nodes, making the distribution more symmetric. The overall pattern is the same in figure 5c except that the initial rapid increase in $[eqn]$ causes a misclassification by our algorithm (see figure 2a).

Second, in figure 5d, $[eqn]$ is clearly successful, whereas $[eqn]$ is less so. As noted above, there is a tendency for $[eqn]$ to reverse its decline near the tipping point, while this tendency is weak. By contrast, as the tipping point is approached, the distribution of $[eqn]$ becomes detectably more symmetric, which successfully increases $[eqn]$ . This change is primarily caused by a faster decrease in $[eqn]$ at the large-degree nodes rather than by that at the small-degree nodes; the latter was the case in figure 5b.

Third, in figure 5f, $[eqn]$ is unsuccessful because $[eqn]$ is bounded in $[eqn]$ for the SIS dynamics. Owing to this restriction, a faster decrease in $[eqn]$ at the large-degree nodes as the tipping point is approached, which occurs in the case of the coupled-well dynamics to let $[eqn]$ perform well (see figure 5d), cannot occur in the case of the SIS dynamics. In this case, the behaviour of $[eqn]$ is apparently not related to that of $[eqn]$ .

Moran’s I

3.3.3.

To explore patterns in the behaviour of Moran’s $[eqn]$ , we show in figure 6 the value of $[eqn]$ , which we refer to as the numerator, $[eqn]$ , which we refer to as the denominator and $[eqn]$ for the same six series of simulations as those shown in figure 5. Note that $[eqn]$ is the ratio of the thus defined numerator to the denominator. The numerator is a covariance function. The denominator is a variance function and the same as $[eqn]$ .

IM and the two quantities defining IM for the simulations shown in figure 5. In each panel, we show the value of IM (orange line), its numerator (blue line) and its denominator (red line). The direction of the simulation sequence is indicated by the arrow beneath the horizontal axis. The classification result for IM is shown in orange text. (a) and (b): coupled double-well, with u as the control parameter. (c) and (d): coupled double-well, with D as the control parameter. (e) and (f): SIS, with D as the control parameter.

In figure 6a,c, both the numerator and the denominator of $[eqn]$ increase in an accelerating manner, which is desirable as an EWS. However, $[eqn]$ , which is the ratio of the two quantities, at first remains around the same average value and then declines as $[eqn]$ further increases. Therefore, the ratio of two apparently successful EWSs generates a poor EWS. In fact, $[eqn]$ is classified as unsuccessful in these two cases. By contrast, in figure 6b,e, the ratio of two increasing, accelerating quantities, which would individually make desirable EWSs, leads to a successful EWS. Conversely, in figure 6d,f, $[eqn]$ is a successful EWS, but neither of its components is. We conclude that the success or failure of $[eqn]$ depends on the intricate balance between the numerator and denominator and that $[eqn]$ ’s performance is not much linked to the performance of the numerator or denominator.

Discussion

We comprehensively analysed the performance of major spatial EWSs across dynamics models, control parameters, simulation directions and networks. We also gave mechanistic insights into the reasons why good performers work well under some simulation conditions but not in others (see figures 5 and 6). Figure 4, which is our main result, indicates that there is no clear overall winner. However, the figure shows some tendency, which we propose as a recommended practice. When the networked complex system shows saddle-node bifurcations (i.e. coupled double-well and mutualistic species dynamics), $[eqn]$ performs the best overall. Furthermore, if the system is known to transit from the lower state to the upper state, $[eqn]$ outperforms $[eqn]$ . However, in ecological and deforestation dynamics, sudden large drops from an upper to lower state, which may reflect a saddle-node bifurcation, are of practical interest, corresponding to mass extinction and deforestation, respectively. Therefore, we recommend $[eqn]$ for these applications. By contrast, when the observables do not show sudden large jumps at tipping points (i.e. SIS and gene-regulatory dynamics), $[eqn]$ is the best at large. However, all the spatial EWSs including $[eqn]$ perform relatively poorly in two of the four simulation conditions without sudden large jumps (i.e. ascending simulations in the SIS dynamics and descending simulations with $[eqn]$ as the control parameter in the gene-regulatory dynamics). Therefore, we conclude that tipping points of these dynamics are relatively difficult to anticipate. We emphasize that $[eqn]$ , which is a commonly used spatial EWS, performs well only under specific simulation conditions (i.e. combination of dynamics with no sudden large jump, $[eqn]$ as the control parameter, descending simulations). Furthermore, computation of $[eqn]$ needs the adjacency matrix of the network, which is not the case for $[eqn]$ , $[eqn]$ or $[eqn]$ .

For the square lattice, some studies report that both spatial correlation and spatial variance anticipate tipping points acceptably well [7,22,27,39,42], whereas others report that spatial correlation outperforms spatial variance [26,37,56] or vice versa [28,43]. Still others report that alternatives such as spatial skewness and spatial kurtosis [23,29], recovery length [18,40,41], information measures [30,46,57] or a combination of measures [16] are more reliable. Other reports say that different spatial EWSs have better or worse performance even given the same dynamics model, depending on, for example, the bifurcation type [20], direction from which the bifurcation is approached [26], or parameter values [13,21,24]. We showed that spatial EWSs show a diversity of results on heterogeneous networks as well. However, we also showed that one of the two best performers on heterogeneous networks (i.e. $[eqn]$ ) performs poorly on the square lattice. Therefore, the findings available for spatial EWSs on the square lattice do not much help us understand how they behave in heterogeneous networks. This article is, to our knowledge, a first comprehensive report for heterogeneous networks. Follow-up studies using other types of spatial EWSs and a larger variety of dynamics models and heterogeneous networks, including examining the effect of network structure, are warranted for future work. Other types of spatial EWSs include spectral reddening [17], mutual information [46], spatial permutation entropy [57] and patch-size distribution and patch shape [7,19].

A side contribution of this article is a new classification method for EWSs into two successful categories and one unsuccessful category. Our proposal was motivated by the insufficiency of the predominantly used performance measure for EWSs, i.e. Kendall’s $[eqn]$ , which was pointed out in the literature [54,55]. Our measure aims to capture how the EWS nonlinearly changes towards a tipping point and regards that an accelerated increase nearer to the tipping point implies a better signal. However, the behaviour of EWSs is probably more diverse than what our classification or Kendall’s $[eqn]$ can capture. We did not provide stopping criteria either, because it is not a focus of the present work. Further work is needed for better assessments of the quality of EWSs, including temporal EWSs.

Our numerical results suggested that well-performing spatial EWSs may depend on the type of tipping points. This possibility is worth further investigation. In real applications, one does not typically know what type of tipping point is being approached. In our opinion, however, a small number of bifurcation types, in particular, saddle-node, transcritical and Hopf bifurcations, seem to be able to represent a large number of empirical systems. Therefore, we propose to focus on investigating these major types of bifurcations to further understand spatial EWSs on networks. The use of various other dynamical systems on networks and the topological normal form of dynamical systems may be fruitful to this end. Furthermore, a recently developed theory that clarified universal features of the interplay between network structure and dynamics on networks (e.g. [58–60]) may be useful for the theoretical investigation of spatial EWSs on networks. Also important is to clarify the behaviour of spatial EWSs in various false positive and false negative scenarios (e.g. when there is a regime shift or discontinuous jumps in $[eqn]$ , but critical slowing down is absent; [54,57,61–63]). We have verified that the present spatial EWSs do not anticipate sudden and large jumps of $[eqn]$ in two major scenarios (described in [54,61]) in which the jumps do not involve critical slowing down. See the electronic supplementary material, S7, for the results. On the other hand, the spatial EWSs can often anticipate continuous bifurcations such as in the SIS dynamics, as expected. Because epidemic outbreaks are certainly a desired application of EWSs [64], the size or continuity of the jump should not be the only main criterion to characterize regime shifts. Furthermore, if a sudden change in the value of a system parameter induces a sudden regime shift (see the electronic supplementary material, S7), by construction, there should be no useful EWS before the parameter value changes. Therefore, we need further discussion on which type of system changes should be successfully anticipated, advancing the existing literature on false positives and negatives of EWSs [54,61].

Although spatial EWSs were originally proposed for spatial regular networks (i.e. the square lattice and its space-continuous limit), the present results suggest that these EWSs are rather more promising for heterogeneous networks, which most complex empirical networks are. Therefore, this work motivates further work on spatial EWSs for complex networks and their applications to empirical data.