STOCHASTIC BIFURCATIONS AND TIPPING PHENOMENA OF INSECT OUTBREAK SYSTEMS DRIVEN BY α -STABLE L´EVY PROCESSES

. In this work, we mainly characterize stochastic bifurcations and tipping phenomena of insect outbreak dynamical systems driven by α -stable L´evy processes. In one-dimensional insect out-break model, we ﬁnd the ﬁxed points representing refuge and outbreak from the bifurcation curves, and carry out a sensitivity analysis with respect to small changes in the model parameters, the stability index and the noise intensity. We perform the numerical simulations of dynamical trajectories using Monte Carlo methods, which contribute to looking at stochastic hysteresis phenomenon. As for two-dimensional insect outbreak system, we identify the global stability properties of ﬁxed points and express the probability density function by the stationary solution of the nonlocal Fokker-Planck equation. Through numerical modelling, the phase portrait makes it possible to detect critical transitions among the stable states. It is then worth describing stochastic bifurcation associated with tipping points induced by noise through a careful examination of the dynamical paths of the insect outbreak system with external forcing. The results give new insight into the sensitivity of ecosystems to realistic environmental changes represented by stochastic perturbations.


Introduction
Ecology takes biology from the relative simplicity of individuals to explain the complexity of interactions between organisms and their environments [7,30].Its significant implications stretch beyond biology into environmental science and the grand challenges facing society [34].The full spectrum of ecological biology encompasses the approaches at the molecular, organismal, population, community and ecosystem levels, as well as relevant parts of the social sciences.
Dynamical system descriptions can help us understand complex natural systems that are only coarsely observed [4,27].The nonlinear dynamics [13,18] of biological phenomena are frequently influenced by unpredictable components due to the complexity and variability of environmental conditions [20].Because Gaussian noise cannot describe intermittent environmental noisy fluctuations, certain complex phenomena are not suitable to be modeled as stochastic differential equations with Gaussian noise.However, a large number of dynamical systems with bursting or intermittent features arise in various scientific fields due to peculiar dynamical features such as the abrupt events in the climate system and the burst-like events in the transcription process of gene regulation.It is more appropriate to model these discontinuous systems associated with a complex structure of the environment using non-Gaussian Lévy processes with càdlàg trajectories.Introducing Lévy noise into the biological insect outbreak system to simulate impulses caused by external disturbances is more near to reality [8,16].Various challenging factors can bring about sudden changes in the number of insects, including many abrupt climatic changes, environmental damages, high-toxic pesticides, etc.The most interesting Lévy noise is given by the α-stable one, arising in local limit theorems for heavy-tailed random walks.For example, coronavirus superspreading is fat-tailed [31].
It is of great practical importance to delve into the underlying dynamics behind in the insect outbreak system as well as other ecological systems [3,11].The insect called spruce budworm is a serious pest and attacks the leaves of the balsam fir tree, which can be a nuisance for some orchardists and farmers.When an outbreak occurs, an impressive amount of harmful budworms can defoliate and kill most of the fir trees in the forest in a few years' time.It is worth mentioning that the budworm population may be affected by sudden environmental noises [22], such as earthquakes, temperature, and hurricanes.We should work harder to discover how habitat loss and climate change are affecting the insects [21,29].
Ludwig et al. [17] analytically modeled the insect outbreak problem by exploiting a slow-fast system: the budworm population evolves on a fast time scale, whereas the trees grow and die on a slow time scale.They qualitatively analyzed the dynamics of a practical model of the interaction between the quickly reproducing budworms and the slowly recovering forest.Grafke and Vanden-Eijnden [12] precisely discussed the slow manifold and bifurcation structures of stochastic insect outbreak model driven by Brownian motion.They numerically investigated the effect of fluctuations on the dynamics, leading to noise-induced transitions between metastable fixed points captured by large deviation theory.
In the current work, we aim to tackle stochastic insect outbreak systems with α-stable Lévy noises.For such complicated systems, stochastic bifurcations and tipping phenomena occur if some external parameters or noise are varied.We would like to explore in detail the rich and subtle interplay of the nonlinear dynamics and the jump behavior in terms of the Lévy noise of the stochastic perturbation.This procedure is a comprehensive treatment both from a technological and scientific point of view, and requires a diversified approach.
This work focuses on the stochastic bifurcations and tipping phenomena for one-and two-dimensional insect outbreak systems affected by α-stable Lévy noises.We partition the bifurcation diagram into different regions "Refuge", "Bistable" and "Outbreak".We especially look at parameter-dependence of the saddle-node bifurcation in one-dimensional insect outbreak model.We revolve around the equilibria, stability and bifurcation of stochastic slow-fast systems.We chiefly analyse nonlinear dynamics and tipping phenomena stemming from the mathematical biology of the tree foliage area and the budworm population in nature.Moreover, we numerically compute system trajectories with Monte Carlo simulations as well as probability density functions based on the nonlocal Fokker-Planck equations in such one-and two-dimensional setups.The results shall be applicable to a wide range of nonlinear ecosystem models under the influence of noise and parameter changes, where stochastic bifurcations and tipping points are of interest.
The remainder of the article is organized as follows.Sections 2 and 3 are the heart of the paper.In Section 2 we consider one-dimensional insect outbreak model given by nonlinear system with α-stable Lévy noise.We also examine saddle-node bifurcation and noise-induced tipping that arise from the key nonlinearity and changing stochastic perturbations.In Section 3 we perform the bifurcation analysis of two-dimensional stochastic insect outbreak system, and describe non-trivial tipping mechanisms in the coupled dynamics of the tree foliage area and the budworm population.Section 4 summarises our findings and presents our conclusions, as well as a number of directions for future study.

One-dimensional insect outbreak model
For a biological system of bifurcation and catastrophe, we now turn to a model for the sudden outbreak of an insect.The proposed stochastic model for the budworm population dynamics is where x ≥ 0 is the spruce budworm, f is the nominal system, L α t is the α-stable Lévy process and σ is the noise intensity.In the one-dimensional insect outbreak system (2.1), f describes a deterministic nonlinear system with of the dynamics has growth rate r, and carrying capacity k that depends on the amount of foliage left on the trees.The predation part x 2 x 2 +1 represents the death rate due to predation, chiefly by birds.The noise term L α t stands for α-stable Lévy process with the index α ∈ (0, 2) of stability, which characterizes the external change of environment.The constant σ ≥ 0 denotes the intensity of noise.An α-stable Lévy motion moves mainly by small jumps if α ∈ (1, 2) is closed to 2, and mainly by big jumps if α ∈ (0, 1] is closed to 0; see [24].When α ∈ (0, 2), the number of the big jumps increases with decreasing α.The large jumps lead to large error terms that cannot be controlled uniformly in time.
In the absence of predators and stochastic fluctuations, i.e., in the case with the help of the method of separation of variables, we calculate the explicit solution The budworm population x(t) grows logistically approaching k.This means the fact that all the individuals are competing for a finite set of resources, and so the growth rate must decrease as the population grows.
When σ = 0, we want to compute the fixed points of the deterministic system ẋ = f (x).Of course, one fixed point always occurs at x * = 0; it is unstable.The intuitive explanation is that the predation is extremely weak for small x, and so the budworm population grows exponentially for x near zero.We are looking for the other nontrivial fixed points of ẋ = f (x), which are given by the solutions of 2) The left-hand side of (2.2) represents a straight line y = r 1 − x k with a x-intercept equal to k and a y-intercept equal to r, and the right-hand side of (2.2) depicts a curve y = x x 2 +1 that is independent of the parameters; their intersections correspond to the nontrivial fixed points for the system ẋ = f (x).The graphs of the line and the curve are shown in Figure 1a.For the fixed parameter r, as we vary the parameter k, the line moves but the curve doesn't.
The key thing to realize is that as we increase k with r = 0.4 fixed, the line rotates counterclockwise about r.However, we can have one, two, or three intersections in (2.2), depending on the value of k.If k is sufficiently small, for instance, k = 5, there is exactly one intersection.While if k is sufficiently large, for instance, k = 15, there exist three intersections a, b and c.Then the fixed points b and c approach each other and eventually coalesce at some critical value around k ≈ 9.561 (in a saddle-node bifurcation) when the line intersects the curve tangentially.After the bifurcation, the line becomes steeper, the only remaining fixed point is a (in addition to x * = 0, of course).Similarly, a and b can collide and annihilate as r is increased.
The deterministic system ẋ = f (x) has at most four equilibria.We recall that x * = 0 is unstable, and also observe that the stability type must alternate as we move along the x-axis.We determine the stability of the nontrivial fixed points: a is stable, b is unstable, and c is stable.Thus, for r and k in the range corresponding to three positive fixed points, the smaller stable fixed point a is called the refuge level of the budworm population, while the larger stable point c is the outbreak level.From the point of view of pest control, one would like to keep the population at a and away from c.The fate of the system is determined by the initial condition x 0 ; an outbreak occurs if and only if x 0 > b.In this sense, the unstable equilibrium b plays the role of a threshold.
An outbreak can also be triggered by a saddle-node bifurcation.If the parameters r and k drift in such a way that the fixed point a disappears, then the population will jump suddenly to the outbreak level c.The situation is made worse by the hysteresis effect, even if the parameters are restored to their values before the outbreak, the population will not drop back to the refuge level.Now we compute the bifurcation curves in (k, r) space where the system undergoes saddle-node bifurcations.But we run into a difficulty since we can not write r explicitly as a function of k (or k explicitly as a function of r).In the calculation, the bifurcation curves will be written in the parametric form (k(x), r(x)), where x determines the position of saddle-node bifurcation point and runs through all positive values.
As discussed earlier, the condition for a saddle-node bifurcation is that the line y = r(1 − x k ) becomes tangent to the curve y = x x 2 +1 .Thus we demand that both the equality (2.2) of the functions and their derivatives (2.4) We substitute this expression (2.4) for r k into (2.2), which allows us to express r solely in terms of x: (2.5) Then inserting (2.5) into (2.4) yields that The condition k > 0 implies that x must be restricted to the range x > 1.This fact can be well understood.Since x represents the tangent point of the curve and the line in Figure 1a, the infimum of x is corresponding to the maximum of the curve where the line is parallel to x-axis, which means that k approaches infinity.Together (2.5) and (2.6) define the bifurcation curves.For each x > 1, we plot the corresponding point (k(x), r(x)) in the (k, r) plane.The resulting curves determine the limiting behavior of r(x) and k(x) as x → 1 and x → ∞.The different regions separated by the bifurcation curves in Figure 1b are labeled according to the stable fixed points that exist.The refuge level a is the only stable state for low r, and the outbreak level c is the only stable state for large r.In the bistable region, both two stable states a and c exist.Characteristic events are near-threshold excitations.When they showup, we should deploy adaptive and prudent strategies, and take more preemptive and precautionary measures to prohibit the situation from getting worse and protect existing trees.Catastrophes are excitations well beyond threshold.We have to deal with the challenges posed by insect outbreak and do what is necessary to prevent wholesale destruction of the forest when the parameters exceed the dangerous threshold value.
By means of Monte Carlo Simulation, several trajectories of stochastic model (2.1) with growth rate r = 0.4 and carrying capacity k = 9.5 are depicted in Figure 2, where the system is monostable with one stable state a.When the noise intensity σ = 0, Figure 2a is effectively commensurate with 50 trajectories of the deterministic system ẋ = f (x).In the presence of noise, from Figure 2b-d with fixed σ = 0.02, the decrease of the index of stability α = 1.99, 1.9 and 1.5, leads to the increase of the relaxation time to the only stable state in comparison to the deterministic case.In other words, the decrease of the index of stability α ∈ (1, 2) results in the increase of the time for the system remaining in the position that a saddle-node bifurcation occurs, and then a large bistable region is produced.
More interestingly still, the orbital relaxation time is shorter as α is smaller and close to 0 in the case α ∈ (0, 1].Only in a very short time the trajectories can be seen if α is too small.This fact was not the same as we had expected in Figure 2a-d.The insect outbreak system is most vulnerable and susceptible to a noiseinduced abrupt change that is nearly instantaneous.When α = 1, the system relaxes to the vicinity of the fixed point with qualitatively different paths before T = 400, as illustrated in Figure 2e.The trajectories might diverge and jump out of the area of concern due to occasional rogue jumps.As seen in Figure 2f with α = 0.5, the flow changes direction.There are many jumps from the bottom to the top before T = 180, which can be seen as the upward tipping from the bottom steady state to the top place x = 60 or other positions, evidenced also by the discrepancy, complexity and unpredictability.The emergence of large divergences makes the original fixed point unattractive such that many paths leave it and escape from one domain of attraction.The impact of Lévy noise shows an interesting stochastic hysteresis phenomenon in insect outbreak system.This is what we refer to as N-induced tipping (N-tipping) associated with a bifurcation, i.e., a loss of stability of a stable state or attractor with respect to changing stochastic perturbations.The clear definition can be found in the Appendix.To be honest, the x = 60 or other positions is not a stable stochastic equilibrium state.However, the bottom steady state near x = 0 becomes unstable in Figure 2f.At this point, the sudden change in stochastic stability indicates the N-tipping.If we fix r = 0.4 and choose the value of k which is great than 9.561, then there are two metastable states of system (2.1) influenced by Lévy noise.The transition from one stochastic equilibrium state to the other one doesn't always happen as we expected.More specifically, several trajectories starting from one of metastable states do not go to the other one as predicted but travel to other places in different ways.In other words, metastable states in stochastic model (2.1) sometimes lost attractiveness rapidly and unpredictably.It conforms to N-tipping again.
Assume that the solution x of stochastic dynamical system (2.1) has a conditional probability density p(x, t|x 0 , 0).For convenience, we drop the initial condition and simply denote it by p(x, t).Using similar calculations as specified in the reference [32], the nonlocal Fokker-Plank equation quantifying the behaviors of system (2.1) is with the initial condition p(x, 0) = δ(x − x 0 ).As this set of nonlocal Fokker-Planck equations is impossible to be solved analytically, they are usually simulated algorithmically; see discussion in the reference [9].We numerically compute probability density function based on the nonlocal Fokker-Planck equation, as shown in Figure 3a.The shape of the probability density does not change much after T = 80.About T = 100, the probability density function approaches the stationary probability density function.The shapes at T = 80 and T = 100 are almost the same, only the peak drops slightly due to the probability churn at the origin.Approximately the unchanged density p(x, t) at T = 100 can be seen as a stable pattern.Figure 3b demonstrates that the deterministic system is monostable, but the probability density function has two peaks, making it analogous to a bistable system.Thus it can be regarded as a bifurcation induced by Lévy noise.
As α = 1.5, 1.1 and 0.8 decreases, the left peak of the stationary probability density function becomes lower and the right one remains the same magnitude.Therefore, the decreasing α leads to the bigger portion of the probability of the right peak.In other words, the insect has a higher probability to be outbreak.We can observe this phenomenon from Figure 4a.Specifically, we examine probability density function in the case α = 0.8, σ = 0.02 and T = 15, 40, 100.The probability flows to the left peak initially and then flows to the right peak, which is clearly presented in Figure 4b.

Two-dimensional insect outbreak system
The two-dimensional stochastic insect outbreak system is obtained by adding α-stable Lévy noise terms as follows: where (X, Y ) ∈ R 2 , and 0 < ε 1 is a small positive parameter measuring slow and fast time scale separation, meaning that in a formal sense where | • | denotes the Euclidean norm.Since dX dt is small except if X 1 − X S = 0, the X-dynamics is called slow in contrast to the fast dynamics in Y .Here the tree foliage area X recovers slowly on time scale O(ε −1 ) to its long time limit S. For the quickly reproducing budworm population Y , the carrying capacity F X of its logistic growth depends on the tree foliage area available.Both populations are subject to stochastic fluctuations.The stochastic terms {L αi t : t ≥ 0}, i = 1, 2 are independent real-valued symmetric α-stable Lévy processes with Lévy triplets (0, 0, ν αi ) on probability spaces (Ω i , F i , P i ), where α i ∈ (0, 2) are the indexes of stability.The positive constants σ i are the intensities of noise.Note that t is a two-dimensional α-stable Lévy process with the Lévy triplet as l = 0 0 , Q = 0 0 0 0 and ν α (du, dv) = ν α1 (du)δ 0 (dv) + ν α2 (dv)δ 0 (du).For the choice ε = 0.05, S = 0.3 and F = 20, a phase portrait of the two-dimensional deterministic system is sketched in Figure 5.The X-nullcline is the line X = 0.3.The Y -nullcline is the geometric shape for which Y 1 − Y 20X − Y 2 X 2 +Y 2 = 0.The nullclines intersect in three fixed points on the plot of two-dimensional vector field showing the position and stability of equilibria.
There exist two distinct negative eigenvalues λ 1 = −0.05 and λ 2 = −0.83,and hence the fixed point Z a is stable.
There is one negative eigenvalue λ 1 = −0.05 and one positive eigenvalue λ 2 = 0.51, and so this fixed point Z b is an unstable saddle point.Evaluating J at the third fixed point Z c = (0.3, 4.7) leads to The eigenvalues λ 1 = −0.05 and λ 2 = −0.57are negative, and then this fixed point Z c is stable.The stable fixed points Z a and Z c provide two different configurations of tree foliage area and budworm density: one with a relatively low presence of budworms, and the other with a budworm outbreak.The insect outbreak model (3.2) exposes bistability between the equilibria Z a and Z c .The slow-fast system (3.2) has a slow manifold M = {(X, h(X)) : X ≥ 0} foliated by the slow variable X.Here h : X → Y is a nonlinear mapping describing the fast variable Y by using the slow variable X as control parameter.In fact, the graph of the nonlinear mapping h is the slow manifold M for the ecosystem (3.2).It means that for a given available tree foliage area X, the budworm population quickly adjusts to a compatible population density h(X).The stable fixed points Z a and Z c lie on the slow manifold M that connects both fixed points.A saddle-node bifurcation occurs at the point Z X = (0.2, 1.9) separating the basin boundary of two attractors, where a stable and an unstable branch of the slow manifold M emerge.
Stochastic bifurcation may be described by a qualitative change of the stationary probability distribution.
Because of nonlocal Fokker-Planck equation for the solution of stochatic model (3.1), the probability density function p(x, y, t) is computed by utilizing the numerical simulations.For a fixed time T = 10, the simulations in Figure 6a-c with the Lévy noise intensities σ 1 = 0 and σ 2 = 0.5 reveal different qualitative behaviors.The decrease of the stability index α 2 = 1.8, 1, and 0.4 leads to sharper peak of top stable state and higher peak of bottom state.This indicates the change of shape of the probability density function p(x, y, t) in response to the stability index α 2 of noise fluctuation on budworms.In the absence of the noise, that is, when σ 2 = 0, the graph of p(x, y, t) is single-peaked as clarified in Figure 6d.If we compare Figure 6c-d, the probability focuses on the vicinity of X−nullcline for σ 1 = 0 and on the vicinity of slow manifold for σ 2 = 0, respectively.Moreover, the transition of the probability follows X−nullcline for σ 1 = 0 or slow manifold for σ 2 = 0, respectively.The Fokker-Planck equation yields quantitative information about the distribution of tipping events.For each t ∈ [0, 100], we plot the time-varying point (X(t), Y (t)) describing the coexistence between the tree foliage area X(t) and the budworm population Y (t) under the influence of noise in the presence of the stability index α 2 = 1.8 with the Lévy noise intensities σ 1 = 0 and σ 2 = 0.05 by the (X, Y ) plane in Figure 7a.We notice that the natural timescale of X(t) is slower than Y (t).It should be pointed out that the noise is only applied to the Y variable accompanied by the noise-free X variable.There exist two possible types of curves with qualitatively different dynamics depending on the settings of X and Y , which are validated by Monte Carlo simulations.To be more specific: there are 4 trajectories of stochastic system (3.1) that display counterintuitive behaviors, where the budworm population Y increases for monotone increasing X but tips to a dangerous and vicious break; and 5 trajectories are exponential decline curves, where the number of Y reduces exponentially due to a slow increase in the tree growth rate of X. From an ecological perspective, the budworms live in the tree foliage area eating leaves [25].If the number of the budworms is small, they can only graze a few leave so that the trees are able to thrive and reproduce.When there are a large number of trees and a small number of budworms, the faster-growing budworms have a sufficient supply of food to eat, and thus the number of budworms increases fast enough.Ultimately, the massive budworms damage and outstrip their food supply.If the increasing insects break out, large quantities of trees will be eaten in a short of time.The plentiful trees is reduced to extremely low numbers.This is an unusual biological phenomenon.The overgrown budworms no longer have enough food to eat, causing them to die of starvation, and their number will decrease amazingly.The populations of budworm and tree can get infinitesimally close to zero and still recover.
Specifically, we fix the stability parameter α 2 = 1.8 and consider non-trivial dynamics of the budworm population Y (t) that varies in time t for the noise strength σ 2 = 0.05 in Figure 7b.Along 4 trajectories undergoing exponential growth, the fluctuating Y (t) becomes noticeably larger and moves to a qualitative collapse of the budworm population past the critical level as t increases.However, the value of Y (t) undergoes a sudden and To further perform a systematic analysis for the effect of the α-stable Lévy process and the strength of noise, we provide an outlook for the interplay between the tree foliage area X(t) and the budworm population Y (t) with the parameters (α 1 , σ 1 , σ 2 ) = (1.8,0.001, 0) in Figure 7c.Noise is not applied to the Y variable.It turns out that a small amount of noise in X cause fluctuations of the 9 trajectories of stochastic system (3.1) with the results of Monte Carlo simulations.There are 5 trajectories that increase too fast, where the explosive growth of the budworm population Y essentially corresponds to insect outbreak.Other 4 trajectories decay exponentially, where Y asymptotically approaches 0.1 after X = 0.25.
Figure 7a together with c effectively showcase that slow-fast system (3.1)appears to be particularly sensitive to the magnitudes of environmental change which are assumed to be α-stable Lévy noises.We detect quantitative and qualitative changes of stochastic system (3.1) by examining the position and stability type of 9 trajectories.Firstly, we see that 4 curves are growing rapidly in regime (α 2 , σ 1 , σ 2 ) = (1.8,0, 0.05), as suggested by Figure 7a.When we consider the stability index α 1 = 1.8 associated with the Lévy noise intensities σ 1 = 0.001 and σ 2 = 0, there exist 5 exponentially increasing curves obtained from Monte Carlo simulations in Figure 7c.These plots illustrate the change of the number of "explosive" trajectories for insect break in different parameter regimes.Secondly, we find that 5 trajectories dwindle at an exponential rate and stay around Y = 0.1 after X reaches 0.23; see Figure 7a.But 4 trajectories exponentially decay to Y = 0.1 and maintain at a low level for a slightly bigger value of X = 0.25, as demonstrated in Figure 7c.
In Figure 7a, noise is not applied to X(t) but stochastic disturbance of the budworm population Y (t) can longitudinally influence the dynamics of the whole model.In Figure 7c, noise is not applied to Y (t) but a small amount of noise in the tree foliage area X(t) cause horizontal fluctuations of the 9 trajectories of stochastic system.For both X(t) and Y (t) subjected to the effects of noises, the paths are more complicated and unpredictable than that of the circumstance where only X(t) (or Y (t)) is catalyzed by stochastic noise.The special case for α = 2 corresponds to a Brownian motion.There are stochastic variations within a narrow range, but the trajectories are continuous without jumps.It is the main difference.
Lastly, there is a strong need to better understand the deterministic dynamics of the budworm population Y .Figure 7d highlights the noise-free Y parameterized by time t ∈ [0, 100] ( the noise intensity σ 2 = 0).We make the following observations: The initial stage of growth of Y (t) along 5 logarithmic curves is approximately exponential, then the growth slows to arithmetic as time goes to t = 25.The value of Y (t) with small variations is in between 4 and 5 after t = 25.Along other 4 paths undergoing exponential decay, the budworm population Y (t) changes from 1.6 to 0.1 through time and holds the line if the time is in close proximity to t = 20.Even more interestingly, external α-stable Lévy noise can drive Y (t) to tipping by comparing Figure 7b and d.The trajectories of Y (t) under the influence of noise present stochastic variations, which can be thought of as Ntipping.Furthermore, the effects of stochasticity accelerate the speed of rising or decreasing budworm population and force Y (t) to move faster, where the value of Y (t) approaches 0.1 at an earlier time t = 10 than expected.
Fortunately, we could go further using similar numerical experiments and probe a myriad of samples to establish a sufficiently accurate result illustrating the effect of noise intensity and the slow manifold on the system behavior through changing of many trajectories.This will lead to the effectiveness and accuracy of science and math simulations in terms of trajectory statistics.

Conclusions and future challenges
This section is dedicated to drawing the conclusion of this work and to presenting future research perspectives.Mathematically, we established the bifurcation analysis and studied N-tipping of one-and two-dimensional insect outbreak dynamical systems in which the noise source was a Lévy process.The environmental changes were suitably assumed to be α-stable symmetric Lévy motions.We conducted numerical simulations of probability density functions and system trajectories, and then formulated the dynamics of the insect outbreak systems.
Let us comment here briefly on possible extensions of those results for further study.N-tipping is intrinsically unpredictable.To predict tipping point is an outstanding and extremely challenging problem.Early-warning signals can be helpful to constrain tipping mechanisms and mathematical biology.It is very meaningful to get early warning signals near some tipping threshold, and the possibility of preventing tipping in real world ecosystems.Quantifying insect outbreak dynamical systems driven by general Lévy processes and examining the effect of external forcing on system's dynamics would be of particular useful for scientific computation and further analysis.In addition, it would be interesting to extend the present additive noise considerations to the case of multiplicative Lévy noise.The viewpoint presented here seems also promising for investigating a more complex biology model, where the jumps multiply the variables.We plan to show those sophisticated contents, simulations and experiments in the future publications.where L α s− is the left limit of L α s at time s.The function ν α (B) = E(N (1, B)) of the Lévy measure is to describe the expected number of jumps of a certain size at time interval (0, 1].Furthermore, define the compensated Poisson measure of L α t via According to the Lévy-Itô decomposition, L α t can be expressed as In particular, L α t − L α s and L α t−s have the same stable distribution S((t − s) 1 α , 0, 0), for t > s ≥ 0. For α ∈ (0, 1), L α t has larger jumps with lower jump probabilities, while for α ∈ (1, 2), it has smaller jumps with higher jump frequencies, and it is a Cauchy process for α = 1.The special case for α = 2 corresponds to a Brownian motion.

A.2 Stochastic bifurcation
At present, there are two main definitions for stochastic bifurcations [1].One is based on the sudden change of shape of the stationary probability density function, i.e., the so-called phenomenological (P) bifurcation; and the other is based on the sudden change of sign of the largest Lyapunov exponent, i.e., the so-called dynamical (D) bifurcation.D bifurcation is a dynamical concept, which is similar to deterministic bifurcations, while P bifurcation is a static concept.Unfortunately, these two definitions do not agree well, and this means that a new definition of stochastic bifurcation need to be explored.
In contrast to bifurcations in deterministic dynamical systems [14], stochastic bifurcations, sometimes associated with noise-induced transitions, are far from well understood.This applies to noises with different statistics and a collective action of different types of external noises.There are few rigorous general theorems and criteria to detect stochastic bifurcations, which are often only verified by computer simulations [19] or for some particular models.While stochastic bifurcations may be characterized numerically by considering the transformation of the shape of the stationary probability distributions of the associated fractional Fokker-Planck equation when the system and noise parameters change, e.g., a transition from unimodal to bimodal (or from bimodal to unimodal) distribution.

A.3 Tipping phenomena
Tipping points are strongly nonlinear phenomena which can be described in layman's terms as large, sudden and often unexpected changes in the state of a system, caused by small and slow changes in the external inputs.The tipping points have been classified, according to whether they involve, predominantly, a bifurcation, noise, or parameter drift.The notion of a tipping point was popularised by Gladwell [10] and has been used in a wide range of applications including ecology [5,26] and climate science [6,15].For example, the tipping point is used by climate scientists to identify vulnerable features of the climate system.If the systems tip, they are likely to have severe impacts on human society.
There are 3 general classes of tipping [28] generated by a family of stochastic non-autonomous systems dx = f (x, µ(t))dt + g(x)dL t , where L represents a Lévy process and µ(t) is a time-varying input.Bifurcation-induced tipping (B-tipping) occurs when system goes through a bifurcation point as µ(t) changes, in which the output from system changes abruptly or qualitatively owing to a bifurcation of a quasi-static attractor.Noise-induced tipping (N-tipping) happens in multistable system when solution switches randomly from one stable state to another, in which noisy fluctuations result in the system departing from a neighbourhood of a quasi-static attractor.Rate-induced tipping (R-tipping) would come into play when system fails to track quasi-static attractor which continuously changes with control parameter µ(t).Furthermore, stochastic nonlinear systems may exhibit tipping phenomena that result from a combination of several of the above [2].

Figure 1 .
Figure 1.(a) Deterministic bifurcation: the blue curve is the graph of y = x x 2 +1 , whereas the purple line represents y = r 1 − x k with r = 0.4, k = 5 in the refuge region; the dashed line shows the saddle-node bifurcation when k ≈ 9.561; the red line with k = 15 depicts bistable behavior.(b) We identify regions with qualitatively different dynamics.The deterministic system ẋ = f (x) is bistable in the middle region.A saddle-node bifurcation occurs on the left boundary of the bistable region.The critical threshold is formed at the right boundary of the bistable region.The outbreak (refuge) region locates on the right (left) of the bistable region.

Figure 5 .
Figure 5. Vector field of the two-dimensional deterministic system: ε = 0.05, S = 0.3 and F = 20.Brown and light green curves denote X-and Y -nullclines, respectively.The three equilibrium points are located where all of the nullclines intersect.