Phase tipping: how cyclic ecosystems respond to contemporary climate.

We identify the phase of a cycle as a new critical factor for tipping points (critical transitions) in cyclic systems subject to time-varying external conditions. As an example, we consider how contemporary climate variability induces tipping from a predator-prey cycle to extinction in two paradigmatic predator-prey models with an Allee effect. Our analysis of these examples uncovers a counterintuitive behaviour, which we call phase tipping or P-tipping, where tipping to extinction occurs only from certain phases of the cycle. To explain this behaviour, we combine global dynamics with set theory and introduce the concept of partial basin instability for attracting limit cycles. This concept provides a general framework to analyse and identify easily testable criteria for the occurrence of phase tipping in externally forced systems, and can be extended to more complicated attractors.

Introduction

Tipping points or critical transitions are fascinating nonlinear phenomena that are known to occur in complex systems subject to changing external conditions or external inputs. They are ubiquitous in nature and, in layman’s terms, can be described as large, sudden and unexpected changes in the state of the system triggered by small or slow changes in the external inputs [1,2]. Owing to potentially catastrophic and irreversible changes associated with tipping points, it is important to identify and understand the underlying dynamical mechanisms that enable such transitions. To do so, it is helpful to consider base states (attractors for fixed external conditions) whose position or stability change as the external conditions vary over time. Recent work on tipping from base states that are stationary (attracting equilibria) has been shown to result from three generic tipping mechanisms [3]:

Bifurcation-induced tipping or B-tipping occurs when the external input passes through a dangerous bifurcation of the base state, at which point the base state disappears or turns unstable, forcing the system to move to a different state [4-6].

Rate-induced tipping or R-tipping occurs when the external input varies too fast, so the system deviates too far from the moving base state and crosses some tipping threshold [7-11], e.g. into the domain of attraction of a different state [12-17]. The special case of delta-kick external input is referred to as shock-tipping or S-tipping [18]. In contrast to B-tipping, R-tipping need not involve any bifurcations of the base state.

Noise-induced tipping or N-tipping occurs when external random fluctuations drive the system too far from the base state and past some tipping threshold [19], e.g. into the domain of attraction of a different state [20-23].

Many complex systems have non-stationary base states, meaning that these systems exhibit regular or irregular self-sustained oscillations for fixed external inputs [24-31]. Such base states open the possibility for other generic tipping mechanisms when the external inputs vary over time. In this paper, we focus on tipping from the next most complicated base state, a periodic state (attracting limit cycle), and identify a new tipping mechanism:

Phase tipping (partial tipping [26]) or P-tipping occurs when a too fast change or random fluctuations in the external input cause the system to tip to a different state, but only from certain phases (or certain parts) of the base state and its neighbourhood. In other words, the system has to be in the right phases to tip, whereas no tipping occurs from other phases.

The concept of P-tipping naturally extends to more complicated quasi-periodic (attracting tori) and chaotic (strange attractors) base states and, in a certain sense, unifies the notions of R-tipping, S-tipping and N-tipping. A simple intuitive picture is that external inputs can trigger the system past some tipping threshold, but only from certain parts of the base state and its neighbourhood. Thus, P-tipping can also be interpreted as partial tipping. Indeed, examples of P-tipping with smoothly changing external inputs include the recently studied ‘partial R-tipping’ from periodic base states [26], and probabilistic tipping from chaotic base states [28,31,32]. Furthermore, P-tipping offers new insight into classical phenomena such as stochastic resonance [20,33,34], where noise-induced transitions between coexisting non-stationary states occur (predominantly) from certain phases of these states and at an optimal noise strength. Other examples of P-tipping due to random fluctuations include ‘state-dependent vulnerability of synchronization’ in complex networks [35], and ‘phase-sensitive excitability’ from periodic states [19], which can be interpreted as partial N-tipping.

Here, we construct a general mathematical framework to analyse irreversible P-tipping from periodic base states. By ‘irreversible’ we mean that the system approaches a different state in the long term. The framework allows us to explain counterintuitive properties, identify the underlying dynamical mechanism, and give easily testable criteria for the occurrence of P-tipping. Furthermore, motivated by growing evidence that tipping points in the Earth system could be more likely than was thought [2,36,37], we show that P-tipping could occur in real ecosystems subject to contemporary climate change. To be more specific, we uncover robust P-tipping from predator–prey oscillations to extinction due to climate-induced decline in prey resources in two paradigmatic predator–prey models with an Allee effect: the Rosenzweig–MacArthur (RMA) model [38] and the May (or Leslie–Gower–May) model [39]. Intuitively, the phase sensitivity of tipping from predator–prey oscillations arises because a given drop in prey resources has distinctively different effects when applied during the phases of the cycle with the fastest growth and the fastest decline of prey. Both the RMA and May models have been used to study predator–prey interactions in a number of natural systems [40-42]. Here, we use realistic parameter values for the Canada lynx and snowshoe hare system [43,44], together with real climate records from various communities in the boreal and deciduous-boreal forest [45].

The nature of predator–prey interactions often leads to regular, high amplitude, multi-annual cycles [46]. Consumer-resource and host–parasitoid interactions are similar, and also often lead to dramatic cycles [47]. In insects, cyclic outbreaks can be a matter of deep economic concern, as the sudden increase in defoliating insects leads to significant crop damage [48]. In the boreal forest, one of the most famous predator–prey cycles is that of the Canada lynx and snowshoe hare [47]. The Canada lynx is endangered in parts of its southern range, and the snowshoe hare is a keystone species in the north, relied upon by almost all of the mammalian and avian predators there [49]. These examples illustrate the ubiquitous nature of cyclic predator–prey interactions, and their significant economic and environmental importance. Their persistence in the presence of climate change is thus a pressing issue.

Anthropogenic and environmental factors are subjecting cyclic predator–prey systems to external forcing which, through climate change, is being altered dramatically in both spatial and time-dependent ways [41,50-54]. In addition to long-term changes due to global warming, there is a growing interest in changes in climate variability on year-to-decade time scales, owing to its more imminent impacts [55]. In particular, increased variability of short-term climatic events manifests itself as, for example, larger hurricanes, hotter heatwaves and more severe floods [53,56-63]. It is unknown how cyclic predator–prey systems will interact with these changes in climate variability.

Beyond ecology, oscillatory predator–prey interactions play an important role in finance and economics [64,65]. Thus, our work may also be relevant for understanding economies in developing countries [66]. Such economies are non-stationary by nature, and it may well be that developing countries have only short phases in their development, or narrow windows of opportunity, during which external investments can induce transitions from poverty to wealth.

This paper is organized as follows. In §2, we introduce the RMA and May models, define phase for the predator–prey oscillations, and describe the random processes used to model climatic variability. In §3, Monte Carlo simulations of the predator–prey models reveal counterintuitive properties of P-tipping and highlight the key differences from B-tipping. In §4, we present a geometric framework for P-tipping and define the concept of partial basin instability for attracting limit cycles. In §5, we produce two-parameter bifurcation diagrams for the autonomous predator–prey frozen systems with fixed-in-time external inputs, identify bistability between predator–prey cycles and extinction, and reveal parameter regions of partial basin instability—these cannot be captured by classical bifurcation analysis but are essential for understanding P-tipping. Finally, we show that partial basin instability explains and gives testable criteria for the occurrence of P-tipping. We summarize our results in §6.

Oscillatory predator–prey models with varying climate

We carry out our study of P-tipping in the context of two paradigmatic predator–prey models, which we present here. We also define ‘phase’ in the context of the predator–prey limit cycles and nearby oscillations. Finally, we introduce our climate variability model.

The RMA and May models

The RMA model [9,38] describes the time evolution of interacting prey and predator populations [67]: In the prey equation, is the low-density (negative) prey growth rate, quantifies the nonlinear modification of the low-density prey growth, the term gives rise to the strong Allee effect that accounts for negative prey growth rate at low prey population density, is the saturation predator kill rate and is the predator kill half-saturation constant. The ratio is often referred to as the carrying capacity of the ecosystem. It is the maximum prey population that can be sustained by the environment in the absence of predators [44]. In the predator equation, represents the prey-to-predator conversion ratio and is the predator mortality rate. Realistic parameter values, estimated from Canada lynx and snowshoe hare data [43,44], can be found in table 1.

Table 1

Realistic parameter values for the RMA model (2.1) and the May model (2.3), estimated from Canada lynx and snowshoe hare data [43,44].

parameter	units	RMA model	May model
r	1/yr	[0,3]	[0,4]
c	ha/(prey⋅yr)	0.19	0.22
α	prey/(pred⋅yr)	800	505
β	prey/ha	1.5	0.3
χ	pred/prey	0.004	n/a
δ	1/yr	2.2	n/a
s	1/yr	n/a	0.85
q	prey/pred	n/a	205
μ	prey/ha	0.03	0.03
ν	prey/ha	0.003	0.003
ϵ	prey/ha	n/a	0.031

As we explain in §2(c), is a piecewise constant function of time that describes the varying climate. This choice makes the non-autonomous system (2.1) piecewise autonomous in the sense that it behaves like an autonomous system over finite time intervals. Therefore, much can be understood about the behaviour of the non-autonomous system (2.1) by looking at the autonomous frozen system with different but fixed-in-time values of .

The RMA frozen system can have at most four stationary states (equilibria), which are derived by setting in (2.1). In addition to the extinction equilibrium , which is stable for , there is a prey-only equilibrium , the Allee equilibrium and the coexistence equilibrium , whose stability depends on and other system parameters: In the above, we include the argument when an equilibrium’s position depends on . The prey and predator densities of the coexistence equilibrium are given by:

The one-parameter bifurcation diagram of the RMA frozen system in figure 1a reveals various bifurcations and bistability, which are discussed in detail in §5(a) and appendix C. Most importantly, as is increased, the coexistence equilibrium undergoes a supercritical Hopf bifurcation , which makes the equilibrium unstable and produces a stable limit cycle . The cycle corresponds to oscillatory coexistence of predator and prey and is the main focus of this study. In the ecological literature, this Hopf bifurcation is referred to as the paradox of enrichment [68]. As is increased even further, disappears in a dangerous heteroclinic bifurcation at , giving rise to a discontinuity in the branch of coexistence attractors. Past , the only attractor is the extinction equilibrium . This heteroclinic bifurcation indicates where complete depletion of the predator becomes part of the cycle. Note that, in the absence of noise, the predator remains extinct once its level reaches zero because the subspace is invariant. Hence the counterintuitive transition to predator extinction at high prey growth rates.

Figure 1

One-parameter bifurcation diagrams with different but fixed-in-time for (a) the autonomous RMA frozen model (2.1) and (b) the autonomous May frozen model (2.3). The other parameter values are given in table 1. (Online version in colour.)

To show that phase tipping is ubiquitous in predator–prey interactions, we also consider another paradigmatic predator–prey model, the May model [39,44]: This model has the same equation for the prey population density as the RMA model, but differs in the equation for the predator population density . Specifically, is the low-density predator growth rate and is introduced to allow prey extinction. In other words, this model assumes that the predator must have access to other prey which allow it to survive at a low density in the absence of the primary prey . The parameter approximates the minimum prey-to-predator biomass ratio that allows predator population growth, and table 1 contains realistic parameter values, estimated from Canada lynx and snowshoe hare data [43,44].

In addition to the extinction equilibrium , which is always stable, the May frozen system has a prey-only equilibrium , an Allee equilibrium , and two coexistence equilibria and , whose stability depends on the system parameters. Further details and analysis of the May frozen model are provided in appendix A.

Phase of the cycle

To depict phase tipping, each point on the limit cycle, as well as in a neighbourhood of the cycle, must be characterized by its unique phase. In the two-dimensional phase space of the autonomous predator–prey frozen systems (2.1) and (2.3), the stable limit cycle makes a simple rotation about the coexistence equilibrium . We take advantage of this fact and assign a unique phase to every point on the limit cycle using a polar coordinate system anchored in : In other words, the phase of the cycle is the angle measured counter-clockwise from the horizontal half line that extends from in the direction of increasing , as is shown in figure 2. Since the values of for the limit cycles in systems (2.1) and (2.3) are three orders of magnitude smaller than the values of , the ensuing distribution of along is highly non-uniform. To address this issue and achieve a uniform distribution of , we include the factor of in (2.4).

Figure 2

Phase portraits showing the (green) predator–prey limit cycles together with their phases and basin boundaries in (a) the autonomous RMA frozen model (2.1) with and (b) the autonomous May frozen model (2.3) with . The other parameter values are given in table 1. Schematic phase portraits depicting all equilibria and invariant manifolds are shown in appendix A, figure 10. (Online version in colour.)

Figure 10.

Schematic phase portraits showing stable (black dots), unstable (black circles) and saddle (black plus signs) equilibria; the stable/unstable manifolds (black dashed curves) of the saddle equilibria; and the (green) limit cycles in (a) the autonomous RMA frozen model (2.1) with and (b) the autonomous May frozen model (2.3) with . The other parameter values are given in table 1. (Online version in colour.)

In the problem of P-tipping, we often encounter oscillatory solutions that have not converged to the limit cycle . Equation (2.4) allows us to define the ‘phase’ of such oscillatory solutions in a neighbourhood of .

Climate variability

Climate variability here refers to changes in the state of the climate occurring on year-to-decade time scales. We model this process by allowing , i.e. the prey birth rate and the carrying capacity of the ecosystem, to vary over time. This variation can be interpreted as climate-induced changes in resource availability or habitat quality. Seasonal modelling studies often assume sinusoidal variation in climate parameters [69-72], but many key climate variables vary much more abruptly [41]. Since our unit of time is years, rather than months, we focus on abrupt changes in climate.

In ecology, abrupt changes in the form of a single-switch between two values of an input parameter are called press disturbances [73].

Guided by the approach proposed in [45,74], we construct a piecewise constant using two random processes; see figure 3a. First, we assume the amplitude of is a random variable with a continuous uniform probability distribution on a closed interval . Second, we assume the number of consecutive years during which the amplitude of remains constant is a random variable with a discrete probability distribution known as the geometric distribution where is a positive integer and . Such an can be viewed as bounded autocorrelated noise. Using actual climate records from four locations in the boreal and deciduous-boreal forest in North America, we choose a realistic value of [45]. We say the years with constant are of high productivity, or Type-H, if their amplitude is greater than the mean . Otherwise we say the years are of low productivity, or Type-L, as indicated in figure 3a.

Figure 3

Results of a Monte Carlo simulation for the RMA model (2.1), where time-varying is generated using and ‘Climate variability’ interval containing . Shown are numerical tipping experiments (B-tipping and P-tipping) for a fixed initial condition . The other parameter values are given in table 1. (a,b) The time profiles of , and in a single tipping experiment. (c) The values of (red) pre and (blue) post each switch that causes a tipping event. (d) States in the phase plane at the time of the switch that causes a tipping event (i.e. states at the ‘tipping time’ defined in definition 4.4), (grey dots) B-tipping and (black dots) P-tipping. (e) The invariant measure of the limit cycle parameterized by the cycle phase . (f ) Probability distribution of tipping phases for (grey) B-tipping and (black) P-tipping. (Online version in colour.)

In the statistical literature, the above form of the geometric distribution models the number of failures in a Bernoulli trail until the first success occurs, where is the probability of success [75].

B-tipping versus P-tipping in oscillatory predator–prey models

In this section, we use the non-autonomous RMA model (2.1) to demonstrate the occurrence of P-tipping in predator–prey interactions. Furthermore, we highlight the counterintuitive properties of P-tipping by a direct comparison with the intuitive and better understood B-tipping.

Note that, in the non-autonomous system, remains the extinction equilibrium, but the predator–prey limit cycle is replaced by (irregular) predator–prey oscillations. Nonetheless, since the system is piecewise autonomous, the dynamics and bifurcations of the autonomous frozen system help us to understand the behaviour of the non-autonomous one.

B-tipping from predator–prey cycles

We begin with a brief description of B-tipping due to the dangerous heteroclinic bifurcation of the attracting predator–prey limit cycle . In the autonomous frozen system, the cycle exists for the values of below , and disappears in a discontinuous way when ; see figure 1a. Thus, we expect one obvious tipping behaviour in the non-autonomous system with a time-varying : To illustrate properties (B1)–(B3), we perform a Monte Carlo simulation of the non-autonomous RMA system (2.1). We restrict the variation of to the closed interval containing the bifurcation point (see the ‘Climate variability’ label in figure 1a, upper arrow), and perform numerical experiments. In each experiment, we start from a fixed initial condition within the basin of attraction of , and let vary randomly as explained in §2c. We allow the system to continue until tipping from predator–prey oscillations to extinction occurs (figure 3b) due to a step change in from to (figure 3a). We then record the values of in red and the values of in blue in figure 3c, the state in the phase space when the switch from to occurs in figure 3d, and the corresponding phase of this state to produce the tipping-phase histograms in figure 3f. B-tipping is identified as the blue dots above in figure 3c, meaning that transitions to extinction occur when changes from to in agreement with (B1) and (B3). The tipping phases corresponding to grey dots in figure 3d, and the ensuing grey histogram in figure 3f, correlate almost perfectly with the green invariant measure of in figure 3e, in agreement with (B2).

B-tipping from predator–prey oscillations to extinction will occur if increases past the dangerous bifurcation level , and the system converges to before switching back to .

B-tipping will occur from all phases of predator–prey oscillations, but phases where the system spends more time are more likely to tip. An invariant measure of can be obtained and normalized to approximate the probability distribution for B-tipping from a phase as shown in figure 3e; see ref. [76] and appendix B for more details on calculating .

B-tipping from predator–prey oscillations cannot occur when decreases over time because does not undergo any dangerous bifurcations upon decreasing .

P-tipping from predator–prey cycles

The most striking result of the simulation is that B-tipping is not the only tipping mechanism at play. It turns out that there are other, unexpected and counterintuitive tipping transitions. These transitions indicate a new tipping mechanism, whose dynamical properties are in stark contrast to B-tipping: Since the unexpected tipping transitions occur only from certain phases of predator–prey oscillations, we refer to this phenomenon as phase tipping or P-tipping.

Tipping from the predator–prey oscillations to extinction occurs when decreases and does not cross any dangerous bifurcations of , which is in contrast to (B1) and (B3). This is evidenced in figure 3c by the blue dots below depicting transitions to extinction when changes from to .

Tipping occurs only from certain phases of predator–prey oscillations, which is in contrast to (B2). This is evidenced by the black dots in figure 3d, and the ensuing black tipping-phase histogram in figure 3f.

The tipping phases do not correlate at all with the invariant measure of shown in figure 3e. This is evidenced by a comparison with the black histogram in figure 3f.

Although P-tipping is less understood than B-tipping, it is ubiquitous and possibly even more relevant for predator–prey interactions. In figure 4, we restrict climate variability in the RMA model (2.1) to a closed interval that does not contain . In other words, we set . Since the time-varying input cannot cross the dangerous heteroclinic bifurcation, all tipping transitions are P-tipping events. Furthermore, owing to the absence of dangerous bifurcations of in the May model (2.3) in figure 1b, P-tipping from predator–prey oscillations to extinction is the only tipping mechanism in figure 5. Note that P-tipping is more likely to occur in the May model, as evidenced by shorter tipping times; compare figures 4c and 5c.

Figure 4

(a,b) and (d) The same as in figure 3 except for taking values from a different ‘Climate variability’ interval that does not contain . As a result, each of the 1000 tipping events is P-tipping. (c) The probability distribution of tipping at time . The other parameter values are given in table 1.

Figure 5.

(a–d) The same as in figure 4 but for the May frozen model (2.3) with taking values from the ‘Climate variability’ interval . Each of the 1000 tipping events is an instance of P-tipping. The other parameter values are given in table 1. (Online version in colour.)

The numerical experiments in figures 4 and 5 serve as motivating examples for the development of a general mathematical framework for P-tipping in §4.

The Allee threshold: intuitive explanation of P-tipping

Intuitively, P-tipping from predator–prey oscillations to extinction in the non-autonomous system can be understood in terms of an Allee threshold in the autonomous frozen system, separating trajectories that lead to extinction from those that approach the predator–prey cycle (figures 2 and 10), and how a given drop in prey resources affects different phases near the predator–prey cycle via the changing Allee threshold.

The shape and position of both the Allee threshold and the predator–prey cycle are modified by a drop in prey resources . The strongest impact is expected when the drop coincides with the region of the fastest decline in prey and a large predator population . These situations occur near the part of the cycle within a range of phases around , which is close to . There, the drop speeds up the prey decline, which, in conjunction with high predation pressure, creates perfect conditions for the ecosystem to move away from the modified cycle, cross the even closer modified Allee threshold and move towards extinction. Indeed, figures 4 and 5 show that P-tipping occurs from a range of phases around . The ecosystem response is very different if the same drop in prey resources coincides with the region of the fastest growth of prey and a small predator population . These situations occur near a different part of the cycle, within a range of phases around , which is away from . There, the drop slows or even reverses the prey growth, but low predation pressure prevents the ecosystem from crossing the distant Allee threshold and helps it adapt to the modified cycle instead. Hence the observed phase sensitivity of tipping from predator–prey oscillations to extinction in the non-autonomous systems.

A geometric framework for P-tipping: partial basin instability

Motivated by the numerical experiments in figures 4 and 5, and the fact that P-tipping is not captured by classical bifurcation theory, the aim of this section is to provide mathematical tools for analysis of P-tipping. Specifically, we develop a simple geometric framework that uses global properties of the autonomous frozen system to study P-tipping from attracting limit cycles and their neighbourhoods in the non-autonomous system. The key concept is basin instability. This concept was first introduced in ([13], Section 5.2) to study irreversible R-tipping from base states that are stationary (attracting equilibria) for fixed-in-time external inputs. Here, we extend this concept to base states that are attracting limit cycles for fixed-in-time external inputs. Our framework will allow us to give easily testable criteria for the occurrence of P-tipping from limit cycles in general, and explain the counterintuitive collapses to extinction in the predator–prey systems from §3.

Not to be confused with the notion of ‘basin stability’ introduced as a measure related to the volume of the basin of attraction [77].

To define basin instability and P-tipping for limit cycles in general terms, we consider an -dimensional non-autonomous system with , and a piecewise constant external input that can be single-switch or multi-switch. When it is important to highlight the dependence of multi-switch inputs on (see equation (2.5)), we write instead of . Note that the RMA (2.1) and May (2.3) models with from §2(c) are examples of (4.1). Furthermore, we write to denote a solution to the non-autonomous system (4.1) at time started from at initial time . We also consider the corresponding autonomous frozen system with different but fixed-in-time values of the external input , and write to denote a solution to the autonomous frozen system at time started from for a fixed .

Ingredients for defining basin instability

One key ingredient of a basin instability definition is the base attractor in the autonomous frozen system, denoted , whose shape and position in the phase space vary with the input parameter(s) . The second key ingredient is the basin of attraction of the base attractor, denoted , whose shape and extent may also vary with the input parameter(s) . For non-stationary attractors , we work with the distance between a solution and the set , and write when this distance tends to zero as . We define as the open set of all whose trajectories converge to forward in time: We often refer to the closure of the basin of attraction of , denoted , which comprises and its boundary, and to the basin boundary of , which is given by the set difference . Additionally, we assume that either all or part of the basin boundary of is a basin boundary of at least two attractors. This property, in turn, requires that the autonomous frozen system is at least bistable, meaning that it has at least one more attractor, other that , for the same values of the input parameter(s) .

The distance between and is .

The third key ingredient is a parameter path , which we define as a connected set of all possible values of the external input . It is important that does not cross any classical autonomous bifurcations of the base attractor .

Definitions of basin instability for limit cycles

In short, basin instability of the base attractor on a parameter path describes the position of the base attractor at some points on the path relative to the position of its basin of attraction at other points on the path. Here, we define this concept rigorously for attracting limit cycles setwise.

Definition 4.1.

Consider a parameter path . Suppose the frozen system has a family of hyperbolic attracting limit cycles that vary -smoothly with . We say is basin unstable on a path if there are two points on the path, , such that the limit cycle is not contained in the basin of attraction of : Furthermore, we distinguish two observable (or typical) cases of basin instability:

We say is partially basin unstable on a path if there are two points on the path, and , such that the limit cycle is not fully contained in the closure of the basin of attraction of , and, for every two points on the path, and , has a non-empty intersection with the basin of attraction of

We say is totally basin unstable on a path if there are (at least) two points on the path, and , such that lies outside the closure of the basin of attraction of

Remark 4.2.

Additionally, there are two indiscernible (or special) cases of basin instability for limit cycles. They cannot be easily distinguished by observation from total basin instability, or from lack of basin instability. However, the indiscernible cases are necessary (although not sufficient) for the onset of partial basin instability and for transitions between partial and total basin instability.

We say is marginally basin unstable on a path if, in addition to (4.2), for every two points on the path, and , the limit cycle is contained in The special case of marginal basin instability separates the typical cases of ‘no basin instability’ and ‘partial basin instability’. Furthermore, it is related to ‘invisible R-tipping’ and to transitions between ‘tracking’ and ‘partial R-tipping’ identified in [26].

We say is almost totally basin unstable on a path if there are (at least) two points on the path, and , such that does not intersect , and, for every two points on the path, and , the limit cycle intersects . The special case of almost total basin instability separates the typical cases of ‘partial basin instability’ and ‘total basin instability’. Furthermore, it is related to transitions between ‘partial R-tipping’ and ‘total R-tipping’ described in [26].

Note that, for equilibrium base states, ‘partial basin instability’ is not defined, whereas ‘marginal basin instability’ and ‘almost total basin instability’ become the same condition.

Guided by the approach proposed in [13], we would like to augment the classical autonomous bifurcation diagrams for the autonomous frozen system with information about (partial) basin instability of the base attractor . The aim is to reveal non-autonomous instabilities that cannot be explained by classical autonomous bifurcations of the frozen system. To illustrate basin instability of in the bifurcation diagram of the autonomous frozen system, we define the region of basin instability of in the space of the input parameters as follows:

Definition 4.3.

In the autonomous frozen system, consider a -smooth family of hyperbolic attracting limit cycles , and denote it with . For a fixed , we define a region of basin instability of as a set of all points in the space of the input parameters , such that is not contained in the basin of attraction of

Partial basin instability and P-tipping

Thus far, we have worked with a loosely defined concept of P-tipping. In this section, we give rigorous definitions of P-tipping for single-switch and multi-switch , show that partial basin instability of for a single-switch is necessary and sufficient for the occurrence of P-tipping from , and discuss the applicability of this result to multi-switch .

Definition 4.4.

Consider a non-autonomous system (4.1) with a piecewise constant input on a parameter path . Suppose the autonomous frozen system has a family of hyperbolic attracting limit cycles that vary -smoothly with .

Suppose is a single-switch that changes from to at time . Suppose also the system is on at . We then say that system (4.1) undergoes irreversible P-tipping from if there are , such that We call a tipping phase associated with each such .

Suppose is multi-switch with a fixed . If leaves the basin of attraction for good, we use to denote the smallest switching time such that We use to denote the corresponding state, and to denote the corresponding tipping phase. We then say that system (4.1) undergoes irreversible P-tipping if, for some initial condition and all realizations of , there are tipping phases and also a non-zero Lebesgue measure subset of that does not contain any tipping phases .

We call the tipping time.

Remark 4.5.

It should be possible to extend definition 4.4 to:

Smoothly varying , for which P-tipping from is expected to depend on the rate of change of [26,30].

Non-periodic attractors such as tori or chaotic attractors, which may require an alternative phase definition. We return to this point in §6.

In general, the occurrence of P-tipping depends on the initial state, the properties of the external input and the topological structure of the phase space. We now show that partial basin instability of for a single-switch is necessary and sufficient for the occurrence of P-tipping from .

Proposition 4.6.

Consider a non-autonomous system (4.1) and a parameter path . Suppose the frozen system has a family of hyperbolic attracting limit cycles that vary -smoothly with , and is partially basin unstable on . Then, for all and , a single-switch parameter change from to gives irreversible P-tipping from if and only if .

Proof.

A single-switch parameter change from to at time reduces the problem to an autonomous initial value problem with initial condition and fixed . It follows from the definition of basin of attraction that only solutions started from are attracted to the limit cycle . Thus, if is partially basin unstable on and , then there will be that give irreversible tipping, and that give no tipping. Conversely, if there is irreversible P-tipping from , then there must be , which implies .

This rigorous statement no longer holds for multi-switch piecewise constant inputs . The reason is that trajectories are no longer guaranteed to converge to the limit cycle , or to the alternative attractor of the frozen system, if the time interval between consecutive switches is short compared with the time of convergence. Additionally, trajectories started in the basin of attraction of may move away from for finite time. These differences allow for two dynamical scenarios that cannot occur in a system that starts on and is subject to a single-switch .

In the first scenario, following a switch, the system leaves the basin of attraction of , but fails to converge to an alternative attractor before the next switch happens, re-enters the basin of attraction of upon the second switch, and avoids P-tipping despite basin instability of . We refer to such events as ‘rescue events’ [45]. Hence, basin instability of for a given switch within a multi-switch does not guarantee the occurrence of tipping upon this particular switch. For the second scenario, we extend the concept of partial basin instability to the whole basin of attraction of . Suppose that is basin stable on , but its basin of attraction is partially basin unstable on . Following a switch, the trajectory moves away from and enters the basin unstable part of the basin of attraction of , then the next switch happens, and the system undergoes P-tipping in the absence of basin instability of . Hence, partial basin instability of need not be necessary for the occurrence of P-tipping with multi-switch .

Keeping in mind that multi-switch P-tipping is defined for all realizations of , it should be possible to show that, for multi-switch piecewise constant :

Partial basin instability of on is sufficient for the occurrence of P-tipping in system (4.1).

If allows trajectories to converge to between all consecutive switches, then partial basin instability of on is necessary and sufficient for the occurrence of P-tipping in system (4.1).

Partial basin instability and P-tipping in predator–prey models

In this section, we start with classical autonomous bifurcation analysis of the predator–prey frozen systems (2.1) and (2.3) to identify parameter regions with bistability between predator–prey cycles and extinction. Then, we show that predator–prey cycles can be partially basin unstable on several parameter paths that lie within these regions of bistability. Finally, we demonstrate that partial basin instability of predator–prey cycles on a path explains the counterintuitive collapses to extinction that occur only from certain phases of predator–prey oscillations, and gives simple testable criteria for the occurrence of P-tipping in the non-autonomous predator–prey system.

Classical bifurcation analysis: limit cycles and bistability

There are four ecologically relevant parameter regions in the predator–prey frozen systems (2.1) and (2.3), shown in figure 6. These regions have qualitatively different dynamics that can be summarized in terms of stable states as follows:

Figure 6.

Two-parameter bifurcation diagrams for (a) the autonomous RMA frozen model (2.1) in the parameter plane, and (b) the autonomous May frozen model (2.3) in the parameter plane. The other parameter values are given in table 1. (Online version in colour.)

Oscillatory Coexistence or Extinction. The system is bistable and can either settle at the extinction equilibrium , or self-oscillate as it converges to the predator–prey limit cycle . Here is where P-tipping may occur; see the green regions in figure 6.

Stationary Coexistence or Extinction. The system is bistable and can settle either at the extinction equilibrium , or at the coexistence equilibrium ; see the yellow regions in figure 6.

Prey Only or Extinction. The system is bistable and can settle either at the extinction equilibrium , or at the prey-only equilibrium ; see the upper pink region in figure 6a.

Extinction. The system is monostable and can only settle at the extinction equilibrium ; see the other pink regions in figure 6.

The region boundaries are obtained via two-parameter bifurcation analysis using the numerical continuation software XPPAUT [78]. This analysis extends our discussion of the one-parameter bifurcation diagrams from figure 1. We refer to appendix C for the details of the bifurcation analysis, and to [79] for more details on classical autonomous bifurcation theory.

Partial basin instability of predator–prey cycles

We now concentrate on the bistable regions labelled ‘Oscillatory Coexistence or Extinction’, apply definitions 4.1 and 4.3 to predator–prey cycles, and show that

Predator–prey cycles can be partially basin unstable on suitably chosen parameter paths.

Both predator–prey models have large parameter regions of partial basin instability. When superimposed onto classical bifurcation diagrams, these regions reveal P-tipping instabilities that cannot be captured by the classical autonomous bifurcation analysis.

Partial basin instability of in the frozen system is sufficient for the occurrence of P-tipping in the non-autonomous system.

The base attractor is the predator–prey limit cycle , and the alternative attractor is the extinction equilibrium . The basin boundary of is the Allee threshold , which can be computed as the stable invariant manifold of the saddle equilibrium : In the RMA frozen model, is the saddle Allee equilibrium , whereas in the May frozen model, is the saddle coexistence equilibrium that lies near the repelling Allee equilibrium . To uncover the full extent of partial basin instability for the predator–prey cycles , we fix a point that lies within the region labelled ‘Oscillatory Coexistence or Extinction’; see figures 7a and 8a. Then, we apply definition (4.7) to identify all points within this region such that the predator–prey limit cycle is not contained in the closure of the basin of attraction of . The ensuing (light grey) regions of partial basin instability bounded by the (dark grey) curves of marginal basin instability are superimposed on the classical bifurcation diagrams in figures 7a and 8a. Note that the basin instability regions depend on the choice of , and are labelled simply for brevity. To illustrate the underlying mechanism in the phase plane, we restrict to parameter paths that are straight horizontal lines from in the direction of decreasing . In other words, we set ; see figures 7a and 8a. When lies on the dark grey curve of marginal basin instability, there is a single point of tangency between and , denoted in figures 7d and 8d. When lies within the light grey region of partial basin instability, there are two points of intersection between and , denoted and in figures 7e and 8e. These two points bound the (red) part of the cycle that is basin unstable. The corresponding basin unstable phases are shown in figures 7b and 8b. Suppose that , and a trajectory of the non-autonomous system is on the same side of as the (red) basin unstable part of . Then, when changes from to , the trajectory finds itself in the basin of attraction of the extinction equilibrium , and will thus approach .

Figure 7

(a) The two-parameter bifurcation diagram for the autonomous RM frozen model (2.1) from figure 6a with the addition of the (grey) region of partial basin instability, for , as defined by (4.7), and the parameter path from . (b) The range of basin unstable phases for the predator–prey limit cycle along . (c–e) Selected phase portraits showing (c) no basin instability for , (d) marginal basin instability for and (e) partial basin instability of on for . Basin stable parts of are shown in green, basin unstable parts of are shown in red. The other parameter values are given in table 1. (Online version in colour.)

Figure 8.

(a) The two-parameter bifurcation diagram for the autonomous May frozen model (2.3) from figure 6b with the addition of the (grey) region of partial basin instability, for , as defined by (4.7), and the parameter path from . (b) The range of basin unstable phases for the predator–prey limit cycle along . (c–e) Selected phase portraits showing (c) no basin instability for , (d) marginal basin instability for and (e) partial basin instability of on for . Basin stable parts of are shown in green, basin unstable parts of are shown in red. The other parameter values are given in table 1. (Online version in colour.)

The striking similarity is that predator–prey cycles from both models exhibit partial basin instability upon decreasing . This decrease corresponds to climate-induced decline in the resources or in the quality of habitat. Furthermore, while the predator–prey cycle in the May model has a noticeably wider range of basin unstable phases, neither cycle appears to be totally basin unstable. All these observations are consistent with the counterintuitive properties (P1)–(P3) of P-tipping identified in the numerical experiments in §3.

Partial basin instability explains P-tipping

Now, we can demonstrate that partial basin instability of in the autonomous predator–prey frozen systems explains and gives simple testable criteria for the occurrence of P-tipping in the non-autonomous systems. The families of attracting predator–prey limit cycles , and their basin boundaries , are the two crucial components of the discussion below.

First, recall the numerical P-tipping experiments from §3, and focus on the crescent-shaped ‘clouds’ of states from which P-tipping occurs; see the black dots in figures 4 and 5. Second, recognize that each P-tipping event occurs for a different value of , and thus from a different predator–prey cycle or its neighbourhood. Therefore, we must consider the union of all cycles from the family along the parameter path bounded by and : which is shown in figure 9. Furthermore, we use the basin boundary of the cycle at the left end of the path to divide into its (light green) basin stable part and (pink) basin unstable part on with . The ‘clouds’ of states from which P-tipping occurs agree perfectly with the basin unstable part of . A few black dots that lie slightly outside the basin unstable part of in figure 9b correspond to those P-tipping events that occur from states that have not converged to the limit cycle and lie visibly away from when the switch that causes tipping happens. Those P-tipping events occur if the time interval during which is shorter than the time of convergence to the limit cycle in the autonomous frozen system. For this particular parameter path, we could not detect any tipping events in the absence of partial basin instability of . However, we could detect multiple ‘rescue events’ described in §4(c) (not shown in the figure). In a ‘rescue event’, the system leaves the basin of attraction of the predator–prey cycle after a switch that gives basin instability, but avoids tipping upon this switch because it re-enters the basin of attraction of the predator–prey cycle after some future switch. ‘Rescue events’ occur if the time interval during which is shorter than the time of convergence to the extinction equilibrium in the autonomous frozen system. In summary, the general concept of partial basin instability of on a parameter path from definition 4.1 is an excellent indicator for the occurrence of P-tipping in the RMA (2.1) and May (2.3) models.

Figure 9.

The concept of partial basin instability on a parameter path with (see definition 4.1) is applied to the union of all predator–prey limit cycles on the path (see equation (5.1)) to explain the counterintuitive P-tipping phenomenon uncovered in figures 4 and 5. The (black dots) states from which the system P-tips to extinction agree perfectly with the (pink) basin unstable parts of for (a) the RMA model (2.1) with , and , and (b) the May model (2.3) with , and . The other parameter values are given in table 1.

Conclusion

This paper studies nonlinear tipping phenomena, or critical transitions, in non-autonomous dynamical systems with time-varying external inputs. In addition to the well-known critical factors for tipping in systems that are stationary in the absence of external inputs, namely bifurcation, rate of change and noise, we identify here the phase of predator–prey limit cycles and nearby oscillations as a new critical factor in systems that are cyclic in the absence of external inputs.

To illustrate the new tipping phenomenon in a realistic setting, we consider two paradigmatic predator–prey models with an Allee effect, namely the RMA model [38] and the May model [39]. We describe temporal changes in the carrying capacity of the ecosystem with real climate variability records from different communities in the boreal and deciduous-boreal forest [45], and use realistic parameter values for the Canada lynx and snowshoe hare system [43,44]. Monte Carlo simulation reveals a robust phenomenon, where a drop in the carrying capacity tips the ecosystem from predator–prey oscillations to extinction. The special and somewhat counterintuitive result is that tipping occurs: (i) without crossing any bifurcations, and (ii) only from certain phases of the oscillations. Thus, we refer to this phenomenon as phase tipping (partial tipping), or simply P-tipping. Intuitively, P-tipping from predator–prey oscillations to extinction arises because a fixed drop in prey resources has distinctively different effects when applied during the phases of the oscillations with the fastest growth and the fastest decline of prey.

Motivated by the outcome of the simulation, we develop an accessible and general mathematical framework to analyse P-tipping and reveal the underlying dynamical mechanism. Specifically, we employ notions from set-valued dynamics to extend the geometric concept of basin instability, introduced in [13] for equilibria, to limit cycles. The main idea is to consider the autonomous frozen system with different but fixed-in-time values of the external input along some parameter path, and examine the position of the limit cycle at some point on the path relative to the position of its basin of attraction at other points on the path. First, we define different types of basin instability for limit cycles, and focus on partial basin instability that does not exist for equilibria. Second, we show that partial basin instability in the autonomous frozen system is necessary and sufficient for the occurrence of P-tipping in the non-autonomous system with a single-switch external input. Furthermore, we discuss the applicability of this result to multi-switch external inputs. Third, we relate our results to those of ref. [26] on rate-induced tipping from limit cycles.

We then apply the general framework to the ecosystem models and explain the counterintuitive transitions from certain phases of predator–prey oscillations to extinction. We use classical autonomous bifurcation analysis to identify parameter regions with bistability between predator–prey cycles and extinction. In this way, we show that predator–prey cycles can be partially basin unstable on typical parameter paths within these bistability regions. Moreover, we superimpose regions of partial basin instability onto classical autonomous bifurcation diagrams to reveal P-tipping instabilities that are robust but cannot be captured by classical bifurcation analysis.

We believe that this approach will enable scientists to uncover P-tipping in many different cyclic systems from applications ranging from natural science and engineering to economics. For example, the predator–prey paradigm is found across biological applications modelling, including epidemiology [80], pest control [81], fisheries [82], cancer [83,84] and agriculture [85,86]. The fundamental relationship described in predator–prey models also appears in many areas outside of the biological sciences, with recent examples including atmospheric sciences [87], economic development [64,65], trade and financial crises [88-90] and land management [91]. External disturbances of different kinds exist in all of these systems, suggesting that the P-tipping behaviours discovered in this paper are of broad practical relevance.

Furthermore, the concept of P-tipping, for base states that are attracting limit cycles with regular basin boundaries, naturally extends to more complicated base states, such as quasi-periodic tori and chaotic attractors, and to irregular (e.g. fractal) basin boundaries [28,31,32,92]. Defining phase for more complicated cycles in higher dimensions, and for non-periodic oscillations, will usually require a different approach. For example, one could define phase for an attracting limit cycle in a multi-dimensional system in terms of its period as a linear function of time . This definition is independent of the coordinate system, and can be extended to every point in the basin of attraction using isochrones. Another approach is to work with a time series of a single observable and use the Hilbert transform to construct the complex-valued analytic signal, and then extract the so-called instantaneous phase [93,94]. This phase variable may provide valuable physical insights into the problem of P-tipping when the polar coordinate approach does not work, or when the base attractor or its basin boundary have complicated geometry and are difficult to visualize. Such systems will likely exhibit even more counterintuitive tipping behaviours, but their analysis requires mathematical techniques beyond the scope of this paper.

Another interesting research question is that of early warning indicators for P-tipping. In the past decade, many studies of noisy real-world time-series records revealed prompt changes in the statistical properties of the data prior to tipping [1,21,22,95], which appear to be generic for tipping from equilibria. However, it is unclear if these statistical early warning indicators appear for P-tipping, or if one needs to identify alternatives such as finite time Lyapunov Exponent [96].