Metastability for discontinuous dynamical systems under Lévy noise: Case study on Amazonian Vegetation

For the tipping elements in the Earth’s climate system, the most important issue to address is how stable is the desirable state against random perturbations. Extreme biotic and climatic events pose severe hazards to tropical rainforests. Their local effects are extremely stochastic and difficult to measure. Moreover, the direction and intensity of the response of forest trees to such perturbations are unknown, especially given the lack of efficient dynamical vegetation models to evaluate forest tree cover changes over time. In this study, we consider randomness in the mathematical modelling of forest trees by incorporating uncertainty through a stochastic differential equation. According to field-based evidence, the interactions between fires and droughts are a more direct mechanism that may describe sudden forest degradation in the south-eastern Amazon. In modeling the Amazonian vegetation system, we include symmetric α-stable Lévy perturbations. We report results of stability analysis of the metastable fertile forest state. We conclude that even a very slight threat to the forest state stability represents L´evy noise with large jumps of low intensity, that can be interpreted as a fire occurring in a non-drought year. During years of severe drought, high-intensity fires significantly accelerate the transition between a forest and savanna state.

crit crit where x is the relative forest cover that grows with the saturating rate G if x > X crit and dies with rate D (assuming G > D > 0). This model has two equilibria: the forest state = − x 1 F D G and the savanna state x S = 0. x F (resp. x S ) exists and is stable if x F > X crit (resp. X crit > 0). X crit is the critical forest cover threshold and directly depends on aridity. Owing to global warming, aridity has a tendency to increase, thus causing a significant displacement of X crit and consequently a size decrease in the basin of the forest state, thus contributing to the x F instability before the perturbations. When X crit reaches the level of x F , the forest state disappears.
In the mathematical literature, this type of dynamical system is known as piecewise-smooth dynamics 18,19 or non-smooth dynamical systems, which are described by differential equations with a discontinuous right-hand side 20,21 . Numerous fundamental questions arise when working with discontinuous dynamical systems. The most basic question is the notion of a solution. Many researchers have contributed to the foundation of this issue 22 , including Filippov, Caratheodory, Krasovskii, Euler and Hermes. However, for our purpose we chose the generalized definition of the solution based on the Filippov theory 21 . Therefore, the AV model can be described by a more general n-dimensional nonlinear system with a discontinuous right-hand side 20 : x ( ) ( , ( )) ( , ( )), , ( , ( )), ,

 
where the initial condition x(0) = x 0 . The right-hand side f(t, x) is assumed to be piecewise-continuous and smooth on  − and  + and discontinuous on the hyper-surface Σ, which is called the switching boundary. The boundary is defined by a scalar switching boundary function h(x). In our model Σ = {x = X crit }; thus, x) is therefore assumed to be C 1 on  ∪ Σ − , and f + (t, x) is assumed to be C 1 on  ∪ Σ + . These functions do not agree at the boundary Σ. The system described by (2) does not define f(t, x(t)) if x(t) is on Σ. One of the methods to overcome this problem is known as Filippov's convex method 21 , which consists of an extension (or convexification) of (2) into the following convex differential inclusion (CDI) (or set-valued extension) F(t, x) is the general case and in our model: , , the convex set with two right-hand side f − and f + , is represented by the following equation: The second fundamental question for this piecewise-smooth dynamical system is the uniqueness of the solution. Obviously, the solution of the initial valued problem (IVP) equation (3) where ∉ Σ x 0 is locally unique because f − (t, x) and f + (t, x) are smooth. Uniqueness problems of IVP may arise when x 0 ∈ Σ or the solution crosses the switching boundary Σ. The solution of the differential inclusion (3) with x 0 ∈ Σ does not satisfy the local uniqueness condition in forward time, as presented in the Methods section. In fact, the type of vector field, defined by (3), is around the switching boundary x = X crit and is called a repulsive (repelling) sliding mode 18 after the solutions diverge from x = X crit , i.e., a solution that starts close to x = X crit will move away from it. However, a solution with the initial condition x 0 = X crit may stay on X crit (because 0 ∈ F(X crit )) or leave the switching boundary by entering either −  or +  . The IVP with the initial condition x 0 = X crit has three possible solutions: Varying the parameters D, G and X crit (i.e. G > D > 0 and x F > X crit ) do not lead the system to a bifurcation. The system experiences structural instability only when X crit reaches the x F level. Thus, considering the objectives of our analysis, we assume that D = 0.2, G = 0.85 and X crit = 0.3. In this case, the convex differential inclusion F(t, x) (3) of the discontinuous AV system (1) is as follows: The discontinuous vector field of the Amazonian vegetation model is shown in Fig. 1(a).

Smooth approximation.
To overcome the problem of non-uniqueness of solutions in the non-smooth dynamical system, and to be able to apply the stability analysis techniques for the metastable states in the case of stochastically perturbed dynamical system, we will achieve a smooth approximation of the discontinuous vector fields. Often, smoothing methods are used to address complicated bifurcations 19 , difficulties in numerical integration 20,24 or problems with the existence and uniqueness of the solution, which arises in differential inclusions with sliding modes. However, this method has some disadvantages, including that it generates stiff differential equations which are numerically expensive to solve, and it does not always lead to the approximations that convey the physical reality 20 in the best way. However, in our case, this approximation is useful, because the smoothness of the vector field is a necessary property in our analysis.
Among the wide variety of smooth approximation methods, the method that stands out is the mollification method. This method was chosen because of the large number of benefits that it provides during the regularization of many ill-posed problems [25][26][27][28][29] , including efficiency, accuracy, robustness, and reduced costs. A detailed explanation of how we carry-out the mollification of (3) is found in the Methods section. Our mollified AV model is now described by the following differential equation with a piecewise-smooth continuous vector field (all computations were performed in Matlab): The mollification of the original AV model changes the dynamical behaviour of the system in the neighbourhood of the discontinuity x = X crit = 0.3, transforming the repulsive sliding mode in the unstable equilibria x A = 0.27 (resp. 0.28, 0.29) for ε = 0.08 (resp. 0.05, 0.025). The vector field f ε (x) and the phase portrait of the mollified AV model are shown in Fig. 1(c,b). 7,11 has provided field-based evidence for the premise that one of the important variables of a landscape, such as tree cover, experiences drastic variations or sharp transitions in responding to climate changes and other stressors. High-intensity fires associated with severe drought may accelerate a widespread degradation of Amazonian forests by abruptly increasing tree mortality 11 . When forest fires do occur under average weather conditions, they typically move slowly, liberating little energy, and they are short in duration and end at night when relative humidity increases. During years of severe drought, the fuel (e.g., leaves, twigs and branches) becomes more abundant and drier, thereby increasing the fires intensity and consequently killing a very high percentage of the trees. An extreme event acts either as an external disturbance, which forest systems can resist, or as a disturbance exceeding the resilience of forest ecosystems and preventing to return to the former dynamical state 8 , thus indicating the presence of the tipping point. Because the interactions between fires and droughts are a more direct mechanism of sudden forest degradation in the south-eastern Amazon 11 , the dynamical model of forest cover evolution goes beyond the typical tipping point, which is modelled via discontinuity in the vector fields; hence, the model must contain a stochastic term that fits perturbations such as fires. Recently, to model these abrupt pulses, burst-like or extreme events have been given higher priority to Lévy perturbations with jumps 12, 13, 30 , because of their properties, such as heavy-tailed distributions, and noncontinuous sample paths. The probabilistic description of the Lévy process is introduced in the Methods section. This description includes a more intuitive but concise premise of the concept that should be understandable to a broad audience. In our stochastically perturbed AV model we incorporate perturbations of only climatic nature, such as climatic changes 31 , aridity 2 , precipitations 3, 4 and fires 11 , which the model may response to by exhibiting metastability. In this way, we consider the fourth case of perturbations, described in the Methods section, which includes the symmetric α-stable Lévy process. The time that the process spends below the value zero is regarded as the hibernation time 32 (as exhibited by plants adapted to a desert environment), in which the trees that are more robust to the aridity can remain in dry state x S = 0 (without dying) for a long time until rain comes.

Randomly perturbed Amazonian vegetation system. Previous research in environmental sciences
Following the above discussion, we include random perturbations of α L t type in system (7) and obtain the following SDE in  1 , where f ε is the mollified drift, and ψ is the noise intensity parameter. According to the existence and uniqueness theorem 12 (Theorem 7.26 p.202) and under the Lipschitz and growth conditions, the SDE (9) has a unique, adapted, cadlag solution X t . The solutions of the deterministic AV system (7) and the stochastically perturbed system (9) for different initial conditions and values of α and ψ are shown in Fig. 2.
The generator A for SDE (9) is For some years, SDEs with discontinuous drift have attracted substantial attention. Different methods and techniques have been proposed to address the difficulties arising from the discontinuities in the vector field. For instance, in ref. 33 the authors have performed a complete qualitative classification for the isolated singular points, i.e., points with deleted neighbourhoods for which the function (1 + |b|)/δ 2 (where b-drift and σ-diffusion coefficient) is not locally integrable. This classification allows for the description of the behaviour of solutions in the neighbourhood of isolated singular points and detects the types of points that disturb the uniqueness of solutions. Stochastic bifurcation analysis has also been studied for different types of models: for smooth and discontinuous oscillators 34 , and for piecewise-smooth ODEs in two dimensions with additive Gaussian noise 35 , among others. The dynamical behaviour of stochastically perturbed solutions near the switching manifold has been studied in refs 36 and 37. Attention has also been paid to the solutions 38,39 , to the transition density function 40 as well as to the noise-induced regularization 41,42 of the SDE with a discontinuous vector field.

Results and Discussion
Stability analysis of the metastable forest state. The stochastic AV model (9) exhibits the metastability phenomenon between the two stable states: the savanna x S and the forest x F . Metastability is defined as a behavioural phenomenon of the solution of the system, which consists of sudden visits with varying durations and frequencies of all domains of attraction 30 . We are particularly interested in performing a stability analysis of the current fertile forest state against stochastic perturbations 2 . To carry out the stability analysis, we study three quantities that provide information on the dynamical behaviour of the system and thus are appropriate for this type of analysis. These are meant to be based on the first exit time, escape probability and stochastic basin of attraction. More information about these quantities can be found in the Methods section.
In the following section, we present the main results of stability analysis for the metastable fertile forest state x F that are based on the Figs 3a-c and 4a-c. By the definitions of the SBA in the case of  1 the basin consists of the two intervals, i.e. thick blue segment and the thick red segment of the well-potential curve see Figs 3c and 4c. This composition directly relates to the two criteria according to which the basin size is defined 14  whose solutions have a "small" probability, measured by the level m set out in Criterion I, of exit from the neighborhood of the fertile forest attractor. In Fig. 3a, that show the escape probability from domain D to D c , the set D I described above remains below the red line (m = 0.5) that is the probability level established by Criterion I.
Going back to Fig. 3c. the thick red segment is the set of initial conditions ∪ = = D D II c i n iII c 1 whose solutions have a "high" probability, measured by the level M set out in Criterion II, of return to the vicinity of fertile forest attractor. In Fig. 3b, that show the escape (return) probability from domain D I c to D I , the set D II c described above remains over the red line (M = 0.7) that is the probability level established by Criterion II. The Fig. 4b,c representing the graphical results of the stability analysis that is based on the second definition of the SBA 43 have the similar description as Fig. 3b,c. The exception is the Fig. 4a that reproduces the mean exit time from D, since in this definition the Criterion I involves mean exit time to define the size of the thick blue segment that is the set of initial conditions D I whose solutions stay longer u(x) ≥ m in the vicinity of the forests attractor. = 1.5, ψ = 1). Given the initial condition x 0 = 0.5, the probability that the forest cover, under the noise influence (α = 1.5, ψ = 1), will undergo the forest-savannah transition is six times higher than the probability of transitions as a result of noise with the parameters α = 0.5 and ψ = 0.1, see Fig. 3a. This result is a clear indication that small noise jumps strongly contribute to x F instability.

The maximum probability value of the forest-savanna transition is detected under Lévy noise with smaller jumps but with higher frequencies and intensities (α
The forest cover remains in the fertile forest state longer under Lévy perturbations with smaller jumps of lower intensities (α = 1.5, ψ = 0.1). As the intensity of the noise increases, the mean residence time of the forest cover in the x F basin decreases, thus contributing to the instability of the forest state. The size of the noise jump inversely influences the mean exit time for various levels of noise intensity: for high noise intensity (say ψ = 1) with the increase of the jump size, which raises the exit time, an enhancement of the forest state stability is seen; however, for a low noise intensity (say ψ = 0.1) with the increase of the jump size, which reduces the exit time, the forest state becomes more unstable, see  basins' definitions are similar. However, these two definitions provide additional information, which is complementary and provides a complete view of the state stability. In the case of the first SBA definition, which is based on the escape probability, the strongest contribution to the size of the basin is due to the set of initial points (0.27, +∞) (resp. (0.32, +∞)) for α = 0.5 (resp. α = 1) whose solutions have the decreased escape probability from the deterministic basin of attraction. However, the entrance probabilities do not significantly contribute to the size of the SBA (see Fig. 3a-c). The smaller basins, i.e. (0.36, +∞) and (0.42, +∞), are obtained for noise with smaller jumps (α = 1, 5) independently of the intensity, and therefore are the cause of the strongest state instability.
In the case of the second SBA definition, on the basis of the mean exit time and escape probability, the time has the same contribution (due to the choice of the criterion m = AMET average mean exit time) to the length of the basin for the different noise parameters α and ψ. What differentiates the length of the basins is the second criterion based on the escape probability, see Fig. 4a-c. Consequently from the second SBA definition, the shorter basins (0.40, +∞) (resp. (0.43, +∞)) are obtained for noise with higher intensity and larger jumps α = 1, ψ = 1 (resp. α = 0.5, ψ = 1). = 0.5, α = 1 and ψ = 0.1). Lévy noise with small jumps (α = 1.5) as well as noise with high intensity (ψ = 1) significantly accelerates the transition between the forest and savanna states, thus causing high instability of the forest. The size of the SBA of x F does not undergo significant variations with the decrease in the values of the mollification parameter ε from the threshold ε = 0.08. This conclusion is confirmed by the results in Figure 5, which contains the SBA of the fertile forest state in the mollified Amazonian vegetation model with different values for the parameter ε = 0.08, 0.05 and 0.025.

Even very small threats to the forest state stability represents Lévy noise with large jumps of low intensity (α
We can explain our results from an ecological point of view, associating Lévy noise with large jumps of low intensity (α = 0.5, 1 and ψ = 0.1) to low-intensity fires that occur in non-drought years. However, during years of severe drought and high-intensity fires, the Lévy noise with small jumps of high intensity (α = 1.5 and ψ = 1) significantly accelerates the transition between the forest and savanna states, thereby causing height instability of the forest.

Conclusion
For the tipping elements in the Earth's climate system the most important issue to address is how stable the desirable state is against random, possibly large perturbations. Thus, we performed a stability analysis of the metastable fertile forest state in a stochastically perturbed Amazonian vegetation model with a discontinuous right-hand side.
Considering that the usual notion of solutions is not suitable for a piecewise-smooth dynamics for our research, we used the generalized definition of the solution from Filippov theory. The AV deterministic system does not define the vector field f(t, x(t)) if x(t) is on a switching boundary (tipping point x = X crit ). To overcome this problem, we extend a discontinuous system into a Fillipov convex differential inclusion. The solution concept, in the sense of Filippov theory, guarantees the existence of a solution for this CDI by the assumption that the set-valued function F(t, x) is upper semi-continuous. The type of vector field, Amazonian vegetation CDI, surrounds the switching boundary x = X crit and is called a repulsive (repelling) sliding mode, for which the uniqueness of solutions is not guaranteed. In fact, the IVP with the initial condition x 0 = X crit has three possible solutions. To overcome the problem of non-unique solutions in the AV non-smooth dynamical system, and to be able to apply stability analysis techniques for the metastable states in the case of stochastically perturbed dynamical systems, we have considered a smooth approximation of the discontinuous vector field. The smooth approximation was been performed by mollification techniques, and we have used the convolution kernel generated by Gaussian function as the mollifier. The mollification of the original AV model changes the dynamical behavior of the system in the neighborhood of the discontinuity x = X crit = 0.3, transforming the repulsive sliding mode into the unstable equilibria x A = 0.27 (resp. 0.28, 0.29) for ε = 0.08 (resp. 0.05, 0.025).
Research in the environmental sciences has provided empirical evidence that tree cover experiences drastic variations or sharp transitions in response to climate changes and other stressors. Recently, to model these abrupt pulses, burst-like or extreme events have been given higher priority to Lévy perturbations with jumps, because their properties, such as heavy-tailed distributions and stochastically continuous sample paths, provide the greatest precision to fit the described phenomena. Therefore, in our stochastic AV model, we have considered perturbations of symmetric α-stable Lévy type, which under the considered conditions, exhibit metastability between the two stable states: savanna x S and forest x F . We have been particularly interested in performing a stability analysis of the current fertile forest state against stochastic perturbations. To perform the stability analysis, we have calculated the following three quantities: mean first exit time, escape probability and stochastic basin of attraction, which provide information about the dynamical behavior of the system and thus are appropriate for this type of analysis.
Our main conclusions include the following three points: the maximum probability value of the forest-savanna transition is detected under Lévy motion with small jumps of high frequency and intensity (α = 1.5, ψ = 1), the forest cover remains in the fertile forest state longer under Lévy perturbations with small jumps of low intensity (α = 1.5, ψ = 0.1); and the fertile forest state is the largest stability basin under the symmetric Lévy noise with the large jumps of low intensity α = 0.5, ψ = 0.1. The results of our analysis also show that even a small threat to forest state stability represents Lévy noise with large jumps of low intensity (α = 0.5, ψ = 0.1). In contrast, a Lévy noise with smaller jumps (α = 1.5) as well as noise with higher intensity (ψ = 1) significantly accelerate the transition between the forest and savanna states, thereby causing high instability of the forest.

Methods
Condition for local uniqueness of the solution of the differential inclusion. The solution of the differential inclusion (3) with x 0 ∈ Σ is locally unique in forward time if (i) The projections of the vector field point to the same side of Σ, i.e. the solution exposing a transversal intersection to the switching boundary: The projections point to Σ, i.e. the solution being an attractive sliding mode: where the normal n(x) perpendicular to a locally smooth switching boundary Σ is given by the gradient of h(x): or if h(x) is non-smooth, by using the generalized differential of h(x) (for more detail see ref. 20, section 2.3): where ∂h(x) is assumed to be bounded.
Mollification method. The general idea of the mollification method is to convolve (i.e. a mathematical operation on two functions that is defined as the integral of the product of these functions after one is reversed and shifted) a discontinuous function with a mollifier (i.e. a smooth function with special properties) to get the piecewise-smooth continuous function, that still remains close to the original generalized function. Mollifier were proposed by Friedrichs 44 in the study of partial differential equations and is defined as follows: A real function η ε is called a Friedrichs' mollifier if 0 with the generator η satisfying the following conditions: Among the possible choices for the generator η we use the Gaussian function η = − π x x ( ) exp( ) Thus, to be able to use the Gaussian kernel as the mollifier, instead of a compact support, it is required that the generator's moments μ k must be finite 28 , i.e., , for all integers (16) k k The convolution operator as well as the mollification are defined as ref. 45: is locally integrable, then its mollification is represented by (17) that is, where B(0,ε) is the closed ball with the center 0 and radius ε > 0,  ⊂ U n is open with the boundary ∂U and The mollified function f ε (x) defined in this way possesses the desired properties (proof can be seen in ref . 45): (i) f ε ∈ C ∞ (U ε ). The mollified function becomes infinitely differentiable; (ii) f ε → f a.e. as ε → 0. The mollified function almost everywhere converges to the original one as the parameter ε shrinks; (iii) If f ∈ C(U), then f ε → f uniformly on the compact subset of U; Lévy perturbations. A brief summary of the main properties of Lévy type stochastic process is useful for the understanding of the analysis method and evaluation of the results. The stochastic process that models the behavior of the landscape variable tree cover has to be a positive process. Thus, we ponder the following fore cases.
i) The First case is the Lévy process absorbed at level 0. Let (X, P x ) be a initial stable Lévy process starting at x > 0 and T = inf{t ≥ 0: X t ≤ 0} is the first hitting time at negative half-line. The probability measure  x is the law under P x of the process  X ( , ) x defined as ≥ .
and it is called killed process or absorbed at level 0. If initial Lévy process (X, P x ) has negative jumps it crosses the level 0 by jumping 46 , so  X ( , ) x vanishes with a jump at 0, i.e. the lifetime of the process is almost surely finite. ii) If (X, P x ) has no negative jumps then  X ( , ) x vanishes continuously at 0, belonging to the Second case of the positive ( + -valued) processes.
iii) The Third case is the Lévy process conditioned to stay positive. The following construction of the law  ⁎ is one of the possible constructions that justifies considering  ⁎ X ( , ) as the Lévy process (X, P x ) conditioned to stay positive: x t x t where t  is the Borel filtration generated by X, i.e.  δ = ≤ X s t ( , ) t s . The infinitesimal generators for these processes were explicitly computed in ref. 46. These three cases of processes can be considered for the AV model in the case of the perturbations such as the pests, diseases 8 , deforestation 2 , etc. For which, in most of cases, the level X = 0 represents absorbing state. iv) The Fourth case is the symmetric α-stable Lévy process. A symmetric α-stable scalar Lévy motion α L t with 0 < α < 2 is defined in a similar way as a Brownian motion, excepting two following properties: i) stationary increments The probability density function for α L t is defined by 1 1 where f α is the probability density function for the standard symmetric α-stable random variable ∼ α X S (1, 0, 0) (for more details see ref. 12).
The generating triplet of α L t is (0, 0, ν α ), with the jump measure, i.e. the expected value of the number of jumps of size dy during the unit time, defined as: 1 where c α is the intensity constant. The jump measure controls the intensity and size of the jumps of the process. So the α-stable Lévy process has finite variation as well as larger jumps with lower jump frequencies for small values of α (0 < α < 1) while it has unbounded variation as well as smaller jumps with higher jump probabilities when α ∈ [1, 2).