## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Predicting tipping points in mutualistic networks through dimension reduction

Edited by Nils Chr. Stenseth, University of Oslo, Oslo, Norway, and approved November 27, 2017 (received for review August 23, 2017)

## Significance

Complex systems in many fields, because of their intrinsic nonlinear dynamics, can exhibit a tipping point (point of no return) at which a total collapse of the system occurs. In ecosystems, environmental deterioration can lead to evolution toward a tipping point. To predict tipping point is an outstanding and extremely challenging problem. Using complex bipartite mutualistic networks, we articulate a dimension reduction strategy and establish its general applicability to predicting tipping points using a large number of empirical networks. Not only can our reduced model serve as a paradigm for understanding the tipping point dynamics in real world ecosystems for safeguarding pollinators, the principle can also be extended to other disciplines to address critical issues, such as resilience and sustainability.

## Abstract

Complex networked systems ranging from ecosystems and the climate to economic, social, and infrastructure systems can exhibit a tipping point (a “point of no return”) at which a total collapse of the system occurs. To understand the dynamical mechanism of a tipping point and to predict its occurrence as a system parameter varies are of uttermost importance, tasks that are hindered by the often extremely high dimensionality of the underlying system. Using complex mutualistic networks in ecology as a prototype class of systems, we carry out a dimension reduction process to arrive at an effective 2D system with the two dynamical variables corresponding to the average pollinator and plant abundances. We show, using 59 empirical mutualistic networks extracted from real data, that our 2D model can accurately predict the occurrence of a tipping point, even in the presence of stochastic disturbances. We also find that, because of the lack of sufficient randomness in the structure of the real networks, weighted averaging is necessary in the dimension reduction process. Our reduced model can serve as a paradigm for understanding and predicting the tipping point dynamics in real world mutualistic networks for safeguarding pollinators, and the general principle can be extended to a broad range of disciplines to address the issues of resilience and sustainability.

Avariety of complex dynamical systems ranging from ecosystems and the climate to economic, social, and infrastructure systems can exhibit a tipping point at which a transition from normal functioning to a catastrophic state occurs (1⇓⇓⇓⇓⇓⇓–8). Examples of such transitions are blackouts in the power grids, emergence of massive jamming in urban traffic systems, the shutdown of the thermohaline circulation in the North Atlantic (9), extinction of species in ecosystems (10⇓⇓–13), and the occasional switches of shallow lakes from clear to turbid waters (14). In fact, the transitions are the consequence of gradual changes in the system, which can, for example, be caused by a slow drift in the external conditions. For example, human activities have caused global warming, leading to a continuous deterioration of the environment and consequently, to species extinction. In an ecological network subject to such habitat changes, nodes and/or links in the network can disappear. As the fraction of disappeared nodes and/or links increases through a critical point, the whole system can reach a point of no return—a tipping point past which the whole network collapses, with all species populations simultaneously becoming zero. To predict the tipping point in complex networked systems is a problem of paramount importance and broad interest.

Given a complex dynamical system, it is of general interest to understand the system dynamics near the tipping point, but often, the high-dimensional nature of the system presents a challenge. In ecological science, a major problem is to determine groupings or appropriate levels of aggregation (15, 16) to enable mathematical and physical analyses of the system to gain insights into the fundamental dynamics while neglecting certain details within the groups or aggregates. For a networked system with a large number of mutually interacting components and many independent parameters, the corresponding phase space dimensionality can be prohibitively high for any direct analysis that aims to gain theoretical insights into the dynamical underpinnings of the tipping point. In such a case, the approach of dimension reduction can turn out to be useful. The purpose of this paper is to apply dimension reduction to a class of bipartite mutualistic networked systems in ecology to arrive at a 2D system that captures the essential mutualistic interactions in the original system. More importantly, it can be used to assess the likelihood of the occurrence of a catastrophic tipping point in the system as the environment continues to deteriorate.

In the development of nonlinear dynamics, dimension reduction has played a fundamental role. For example, the classic Lorenz system (17), a system described by three ordinary differential equations (ODEs) with a simple kind of nonlinearity, is the result of drastic reduction in dimension from the Rayleigh–Bénard convection equations with an infinite phase space dimension. Study of the reduced model can lead to insights into dynamical phenomena not only in the original system but also, beyond. In this sense, the reduced model may be said to possess certain features of universality. With regard to tipping point dynamics in complex networked systems, a recent work (18) treated mutualistic, bipartite networked systems and derived a 1D reduced model. In particular, pollinator–plant networks in nature can be regarded effectively as a bipartite network, where any direct interaction in the network is between a pollinator and a plant. Pollinators or plants among themselves, of course, are also connected, albeit indirectly, where two pollinators are regarded as connected if they interact with the same plant, and the same applies to the plant–plant connections. This results in two projection networks: one for the pollinators and another for the plants. The reduced 1D model (18) applies then to either of the projection networks. It was shown that the 1D model can lead to a resilience function, a function describing the emergence of a tipping point as a properly normalized system parameter is changed continuously. The resilience function was speculated to be universal in that it resembles the actual functions obtained from a number of empirical pollinator–plant networks.

While the 1D model is simple and amenable to analysis, it is from a projection network either of the pollinators or of the plants. In the dimension reduction process, certain features of the most fundamental dynamical property of the original bipartite network are lost: mutualistic interactions. To take into account these interactions, a reduced model needs to be simple but not simpler: a 2D model is necessary to capture the bipartite and mutualistic nature, with one collective variable for the pollinators and another for the plants. Consequently, as a single parameter of the system is varied, one can define two resilience functions: one for the pollinators and another for the plants.

We proceed through a series of steps to analyze the data and predict tipping points. (*i*) We develop a method to obtain, for mutualistic networks of arbitrary size, a set of two nonlinear ODEs, where the dynamical variables are the average abundances of the pollinators and plants, respectively. The 2D system contains two key parameters that can be fixed for any given real bipartite mutualistic network. (*ii*) From the generic 2D model, we calculate the average abundances both of the pollinators and of the plants, each as a function of two key parameters that capture the variations among different empirical networks. Each function gives a 2D surface in the 3D space of the average abundance and the two parameters. Since for each real network, the values of the two parameters are fixed, we calculate the values of the average abundances. The remarkable finding is that, for the 59 available plant–pollinator networks from real data recorded over the world, all of the actual average abundance values fall on the 2D surface obtained from the reduced model, providing support for its validity and generality. (*iii*) We calculate, for each real network, two resilience functions by considering two types of parameter variations: (*a*) the fraction of removed pollinators and the associated fraction of removed links and (*b*) the decay rate of individual species. The motivation behind the choice of these parameters is that, because of the continuous deterioration in the environment as a result of, for example, human activities, certain pollinators would disappear and so would the associated mutualistic interactions. In the case where a species manages to survive, the increasingly hostile environment makes it difficult to be sustained, leading to an increase in its decay rate.

The approach finally leads to a method for predicting tipping points. For each real network, we compare the resilience functions with those from the corresponding reduced 2D model and find a good agreement (even in the presence of stochastic disturbances), indicating that the 2D reduced model captures the essential dynamics of the real systems and can thus be used for probabilistic prediction of the occurrence of a tipping point as some parameters reflecting the environmental deterioration change.

## Results

### Nonlinear Dynamical Networks of Mutualistic Interactions and a General Reduced Model in Two Dimensions.

We investigate all mutualistic pollinator–plant networks available from the Web of Life database (www.Web-of-Life.es). There are altogether 59 networks, which cover a wide geographic range across different continents and climatic zones. The structures of the networks are quite different from each other, as is the number of species in each network. Despite the differences, the network dynamics can be described by a set of first-order, nonlinear ODEs, with the total number of equations (the phase space dimension) being the number of species in the network (both pollinators and plants) (13, 19). Considering a generic nonlinear dynamical system described by such ODEs with arbitrary numbers of pollinators and plants, we articulate a dimension reduction process to obtain an average system in two dimensions. As will be shown, the 2D system can capture the essential dynamical features of all 59 real mutualistic networks.

A generic mathematical model for mutualistic interactions includes the following processes: intrinsic growth and intraspecific and interspecific competition as well as mutualistic effects of plants and pollinators. We use the letters P and A to denote plants and pollinators, respectively. Let *i*th plant and the *i*th pollinator, respectively; α is the intrinsic growth rate in the absence of intraspecific and interspecific competition as well as any mutualistic effects. The factors that affect the intraspecific and interspecific competition, such as light and nutrients for the plants at the breeding sites for animals, are characterized by the parameters

In *SI Appendix*, *Note 1*, we detail the steps of our dimension reduction procedure, which leads to the reduced model*SI Appendix*, *Note 1*): (*i*) unweighted average, (*ii*) degree-weighted average, and (*iii*) eigenvector-based average.

Does our reduced 2D system capture some basic properties of real bipartite mutualistic networks? Treating the effective mutualistic interaction strengths (*A* for the pollinator and in Fig. 1*B* for the plant. We can then calculate, for each real network, the parameters *SI Appendix*, Table S1) are near the corresponding surfaces from the reduced model, providing preliminary evidence that the reduced model captures the essential behavior of the real networks from a wide geographical range across continents and climatic zones. As we will show, however, the detailed averaging process can play a role in the model’s predictive power of the average abundances and the tipping point.

Reducing the high-dimensional system (*i*) to the effective 2D system (*iii*) entails an inevitable loss of detailed information about the original system. However, since predicting the occurrence of the tipping point through the system resilience function is our goal, the primary question is whether the reduced 2D system has predictive power, despite the loss of certain details about the dynamical evolution of the original system. In the following, we present strong evidence that the answer to this question is affirmative.

A resilience function is a relationship between the average species abundance and some parameter with variations that reflect the impact on the environment caused by, for example, global warming or direct human activities, such as overuse of pesticides (21⇓–23), where a larger impact corresponds to a higher value of the parameter. Since pollinators are more vulnerable to environmental changes than plants, we focus on two parameters: (*i*) *ii*) κ—the average pollinator decay rate. From the standpoint of plants, the disappearance of a specific pollinator means the loss of a number of links, as any pollinator typically interacts with several plants. Thus, we will also consider the parameter *i*) characterizes the pollinator decay caused by a decrease in the pollinator growth rate and/or an increase in the pollinator mortality rate. Continuous deterioration of the environment leads to a gradual increase in the average decay rate κ. While we have studied all 59 mutualistic networks derived from real data, we report the detailed validation results from 2 representative networks: network A obtained from data recorded at Tenerife, Canary Islands (24) and network B from an empirical study at Hestehaven, Denmark (25). Network A has 38 pollinators and 11 plants, and there are 106 mutualistic interactions. Network B has 42 pollinators and eight plants, with 79 mutualistic interactions. The structures of the two networks are shown in Fig. 2 *A* and *C*, respectively, whereas their matrix representations are shown in Fig. 2 *B* and *D*, respectively.

### Resilience Functions in Systems Without a Tipping Point.

We first examine the case where the system does not exhibit any tipping point. For each value of *SI Appendix*, Fig. S1.) Fig. 3*A* shows, for network A, four types of pollinator abundances vs. *A*, red) and three from the reduced model with different averaging methods (Fig. 3*A*, blue, black, and cyan corresponding to averaging methods *i–iii*, respectively). We see that averaging method *i* leads to abundance variations that are in good agreement with those from the original network, while systematic deviations exist for the results from averaging methods *ii* and *iii*, although they agree with each other. Fig. 3*B* shows, for network A, the corresponding plant abundance vs. *B*, green. The averaging process *i* leads again to average abundance variations in agreement with those from the original network. The results for network B are shown in Fig. 3 *C* and *D*. We see that, for both networks in the parameter setting studied, the remaining pollinators and plants never come close to extinction, even when the fraction of removed pollinators approaches one. The reason is precisely mutualistic interactions: even if there is only one remaining pollinator, at least one plant will be connected with this pollinator. Because of the mutualistic interactions, both the pollinator and plant will survive, and the network system does not exhibit a tipping point. In this parameter regime, the small migration rate has no effect on system dynamical features, such as the absence of a tipping point. The main message of this example is that the reduced model is capable of capturing the abundance variation patterns of both pollinators and plants in the original networked systems.

A phenomenon in Fig. 3 is that, as

### Power of the Reduced Model in Predicting Tipping Points.

We now consider parameter regimes where the mutualistic network system exhibits a tipping point. An examination of the individual resilience functions (*SI Appendix*, Fig. S2) for networks A and B reveals that, as *i*, the abundance values are somewhat smaller than those from the original system, while the opposite behavior occurs for the reduced model with averaging method *ii* or *iii*. These deviations are expected, considering the drastic approximations used in deriving the reduced model. While all three types of average in the reduced model generate results indicating the occurrence of a tipping point, a key issue is the accuracy of the predicted parameter value where a tipping point is reached. In particular, if a reduced model does indeed possesses predictive power for the tipping point, it should predict its occurrence at the correct critical point in the original network. In this regard, we see that the reduced model with averaging method *i* fails to predict the location of the tipping point, whereas the 2D model with averaging method *ii* or *iii* yields the true critical point. Since methods *ii* and *iii* are based on some sort of weighted averaging process (i.e., with respect to nodal degrees and eigenvectors), the results in Fig. 4 point at the importance of imposing weighted averaging process in the dimension reduction process. This conclusion holds not only for the two networks in Fig. 4 but also, for other networks studied.

While structural change in the network, such as gradual removal of nodes (pollinators), can trigger a tipping point, there are alternative scenarios, such as parameter changes. Here, we consider the situation where the mutualistic network structure remains intact, but the death rate of the pollinator increases because of environmental deterioration. Specifically, we increase the pollinator decay rate κ from zero and calculate the species abundances from the original network and from the three averages of the reduced system. The results are shown in Fig. 5. We see that, similar to the case where the network structure is altered through continuous removal of pollinators (compare with Fig. 4), the reduced model through averaging method *ii* or *iii* has a remarkable predictive power for the tipping point in that the predicted critical value of κ at which the species abundances collapse to zero agrees well with that from the original system.

In our computations, we set a relatively small value for the migration rate for all pollinator species:

### Role of Network Randomness in the Reduced Model.

Our extensive computations of the large number of empirical mutualistic networks indicate that a 2D reduced system obtained through degree- or eigenvector-weighted average can correctly predict the tipping point, while the reduced system with unweighted averaging fails to do so. One plausible reason is the lack of sufficient randomness in the network structure. In fact, despite the large variations in their size and structure, the real networks are not quite as random. For a purely random network, either unweighted or weighted averaging has the same effect on the reduced model. To test this proposition and to further show the importance of weighted averaging in the dimension reduction process for real networks, we study artificial random mutualistic networks. Concretely, we construct an ensemble of model networks of 80 pollinators and 40 plants with 960 mutualistic links randomly distributed between the pollinators and the plants and obtain a 2D system through the three averaging processes. As shown in Fig. 6, not only are the three averaging processes able to predict accurately the tipping point, but the differences between the effective abundances that they produced and the actual abundance are small.

### Robust Predictive Power of the Reduced Model Against Stochastic Disturbances (Noises).

Is our reduced model capable of predicting the tipping point when stochastic disturbances are present in the original network? To address this question, we test the predictive power of the reduced model by considering stochastic abundance and parameter fluctuations. First, we assume that there is additive, independent white Gaussian noise in the dynamic equation for the abundance of each species. The results are shown in Fig. 7, where the color legends are the same as those in Fig. 3. We see that, despite the additive noise, the reduced model still predicts correctly the tipping point, where the performance of the model derived using the degree- and eigenvalue-averaging methods is better than that of the model based on unweighted averaging. Second, we study the case where there is randomness in the intraspecific competition rate as motivated by the consideration that, in reality, the intensity of intraspecific competition varies from one species to another. The results are shown in Fig. 8. We see that, despite the large variations in the species abundances caused by the parameter perturbation, the reduced model based on a weighted average method (averaging method *ii* or *iii*) is still capable of predicting the correct tipping point as in Fig. 4.

### Effects of Interspecific Competition.

So far, we have neglected interspecific competitions, as they are generally much weaker than intraspecific competition. Mathematically, interspecific interactions can be modeled through nonzero off-diagonal elements in the competition matrices **1**, which change the structures of these matrices. It is thus useful to investigate the effects of interspecific competition. Our computations reveal that the 2D reduced model captures all essential features of the mutualistic networked systems, even in the presence of interspecific competition, as shown in Figs. 9 and 10, for the cases where tipping points are absent and present, respectively.

Fig. 9 is obtained under the same setting as that of Fig. 3, except that interspecific competition is now included. Comparing Fig. 9 with Fig. 3, we see that the competition results in somewhat lower pollinator and plant abundance, which is intuitive. For small values of

Fig. 10 shows the effects of interspecific competition in the parameter regime where there is a tipping point. Comparing Fig. 10 with Fig. 4, we see that, with the inclusion of interspecific competition, the 2D reduced model derived from the degree- or eigenvector-weighting scheme predicts the tipping point accurately, which is similar to the case where such competition is absent. The mere effect of a reasonable amount of interspecific competition is simply reduced abundances. (For sufficiently large values of the interspecific interaction strength, both the original and the reduced 2D models agreeably give the result of species extinction.)

### Occurrence of a Tipping Point in the 2D Parameter Space.

Both structural change (removal of certain pollinator species) and parameter change (increase in κ) can lead to a tipping point. It is useful to examine the occurrence of a tipping point in the 2D parameter space *SI Appendix*, Fig. S3 for

### A Mathematical Analysis of the System Dynamics with an Explanation of the Emergence of a Tipping Point.

The 2D reduced model provides mathematical insights into the emergence of a tipping point. The relevant quantities are the steady-state values of the species abundances. Setting *SI Appendix*, *Note 3*).

We first consider the parameter regime in which the system does not exhibit a tipping point (e.g., Figs. 1 and 3). In this case, we have **4** can be conveniently expressed in terms of the following algebraic equation for **5**, one is positive, and another is negative. The abundance values of

We next consider the parameter regime, in which the mutualistic system exhibits a tipping point (Fig. 4) (i.e., *SI Appendix*, *Note 3* for an initial state with high abundances, we have *SI Appendix*, Eqs. **S14** and **S16**. After the tipping point has been reached, the steady-state solutions are given by *SI Appendix*, Eqs. **S15** and **S17**. In this case, the physically meaningful solutions are

For a tipping point induced by an increase in the species decay rate κ (Fig. 5), we have **5**. For larger values of κ, the solutions of Eq. **5** become complex, which are physically unrealistic. The condition for the onset of complex solutions is

As the value of κ approaches *SI Appendix*, Eqs. **S14** and **S17**. In particular, we have

Overall, the 2D reduced model thus provides a mathematical paradigm by which a number of distinct dynamical phenomena in mutualistic interacting networks can be understood. For example, as pollinators are removed one after another from the system (i.e., as

The analysis leads to insights into the effect of varying the value of the parameter h, the half-saturation constant. The reason that we choose different values of h for different parameter setting is to ensure that the system exhibits a tipping point as the value of κ or f is varied. From the mathematical analysis, for cases where a tipping point occurs, h will affect the critical values of κ,

## Discussion

Complex dynamical systems exhibiting a tipping point are widespread, and it is of interest to understand the dynamical mechanism of the tipping point and to develop predictive tools. To accomplish these goals, a viable solution is dimension reduction. We focus in this paper on bipartite mutualistic networks, not only as a concrete example to show the use of dimension reduction but also, because of the fundamental values of safeguarding pollinators to human survivability (28). In a mutualistic network system, a tipping point typically exists. As the environment continues to deteriorate, the system can drift toward the tipping point, where the catastrophic phenomenon of pollinator collapse will occur. The backbone that supports the functioning of such a network is mutualistic interactions between the pollinators and plants. To understand the role of the interactions with respect to the emergence of a tipping point, both species of the bipartite network must be retained in a reduced model. That is, the minimum dimension of the reduced system should be two [a 1D reduced model (18) is inadequate to describe mutualistic interactions]. With this in mind, we carry out a dimension reduction process by resorting to different types of averaging methods for species abundances. In particular, given an empirical mutualistic network, we carry out averaging processes to arrive at a 2D model with two collective dynamical variables: one for the pollinators and another for the plants. The average can be either unweighted or weighted. We show that our 2D reduced model captures the essential features of all 59 available real world mutualistic networks, not only in terms of the average abundances but more importantly, in terms of the occurrence of the tipping point, even in the presence of stochastic disturbances. We also find that, because of the lack of sufficient randomness in real mutualistic networks, a weighted average (e.g., based on degrees or eigenvectors) is necessary for the reduced model to exhibit a tipping point at the same critical parameter value as with the original network. Our 2D model can thus serve as a generic paradigm for understanding the tipping point dynamics in real world mutualistic networks. For example, the 2D model can be exploited to investigate a variety of nonlinear dynamical phenomena in mutualistically interacting networked systems, such as bifurcations (29), basin structures (30, 31), crises (32), and transient chaos (33⇓⇓⇓⇓–38), which would otherwise be infeasible with the original systems because of their high dimensionality.

## Acknowledgments

This work was supported by National Science Foundation Grant 1441352. Y.-C.L. acknowledges support from the Vannevar Bush Faculty Fellowship Program sponsored by the Basic Research Office of the Assistant Secretary of Defense for Research and Engineering and funded by Office of Naval Research Grant N00014-16-1-2828 W.L. is supported by the National Science Foundation of China under Grant 61773125.

## Footnotes

- ↵
^{1}To whom correspondence may be addressed. Email: amhastings{at}ucdavis.edu or Ying-Cheng.Lai{at}asu.edu.

Author contributions: J.J., A.H., and Y.-C.L. designed research; J.J. performed research; J.J., Z.-G.H., T.P.S., W.L., C.G., A.H., and Y.-C.L. analyzed data; and Y.-C.L. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

See Commentary on page 635.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1714958115/-/DCSupplemental.

- Copyright © 2018 the Author(s). Published by PNAS.

This open access article is distributed under Creative Commons Attribution-NonCommercial-NoDerivatives License 4.0 (CC BY-NC-ND).

## References

- ↵
- Gladwell M

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Lontzek TS,
- Cai YY,
- Judd KL,
- Lenton TM

- ↵
- Gualdia S,
- Tarziaa M,
- Zamponic F,
- Bouchaudd JP

- ↵
- ↵
- ↵
- Dai L,
- Vorselen D,
- Korolev KS,
- Gore J

- ↵
- Tylianakis JM,
- Coux C

- ↵
- ↵
- Scheffer M

- ↵
- ↵
- ↵
- ↵
- ↵
- Rohr RP,
- Saavedra S,
- Bascompte J

- ↵
- Dakos V,
- Bascompte J

- ↵
- Goulson D,
- Nicholls E,
- Botías C,
- Rotheray EL

- ↵
- Kerr JT, et al.

- ↵
- Koh I, et al.

- ↵
- ↵
- Montero AC

- ↵
- ↵
- ↵
- ↵
- Ott E

- ↵
- ↵
- ↵
- ↵
- Hastings A,
- Higgins K

- ↵
- ↵
- ↵
- ↵
- Hastings A

- ↵
- Lai YC,
- Tél T

## Citation Manager Formats

## Sign up for Article Alerts

## Article Classifications

- Biological Sciences
- Ecology

- Physical Sciences
- Sustainability Science

## See related content: