Many ecosystems are characterized by a rich nonlinear behaviour including multistability as well the appearance of sustained oscillations and catastrophic shifts [Scheffer_2001, Rietkerk_2004, Levin_1998, Jorgensen_2007, Capit_n_2010]. For example, catastrophic/regime shifts have been observed in shallow lakes[Scheffer_2007], the North Pacific, in the North Sea and Caribbean coral reefs [deYoung2008], in regions of northern Africa going abruptly from vegetated to desert conditions.
Frequently, the disturbances in ecosystems that drive such regime shifts are caused by human interactions [Ling_2009, Genkai_Kato_2006, BURNEY_2005, Gordon_2008]. Thus, the understanding from a mathematical point of view of how disturbances/and or state changes influence the dynamics of ecosystems is very important in order to prevent and control ecological disasters [Scheffer_2001, Villa_Mart_n_2015]. From a nonlinear dynamics perspective, a sudden shift from one regime to another can be observed if the system has different coexisting stable regimes [Scheffer_2001, Innes_2013, Villa_Mart_n_2015, Henderson_2016, Russo_2019]
. In such cases, as a consequence of a vector state perturbation, the system may jump to a different stable state due to the crossing of its basin of attraction boundary[Villa_Mart_n_2015]. Another way, in which sudden transitions may be observed is by parameter perturbations beyond a bifurcation point [Villa_Mart_n_2015]. Thus, many studies in the field pinpoint the significance of bifurcation theory to explain, analyse and ultimately forecast beforehand sudden regime shifts [Kooi2003, Scheffer_2001, Scheffer_2010, Bauch_2016, K_fi_2012, Ling_2009, Henderson_2016, Engler_2017, Cimatoribus_2013, Dijkstra_2009, Troost_2007]. However, despite the well-established tools of bifurcation theory, many studies still use simple temporal simulations and/or simple linear numerical analysis to study models that are used to approximate real-world ecosystems. Recently, we have stressed the importance of using the arsenal of analytical and numerical methods of bifurcation theory and particularly that of codimension-2 bifurcation analysis to systematically identify and characterize criticalities in ecological models and in general in physicalm biological and natural systems [Russo_2019]. However, only few studies have used both analytical and numerical bifurcation analysis (i.e. not just temporal numerical simulations for the verification of the analytical findings, which do not suffice to identify and trace in the parameter space branches of far-from-equilibrium bifurcation points such as limit points of limit cycles), thus to study in a combined and complete way not only bifurcations at equilibrium but also bifurcations far-from-equilibrium (see for example [Fujii_1982, Kooi2003, Troost_2007, Hamzah_2007, Dai_2010, Yu_2014]).
Here, we provide both analytical and numerical bifurcation results for a forest-grassland model with human interaction proposed by Innes et al. [Innes_2013]. Forest-grassland mosaic ecosystems are typical examples of ecosystems where two species (forest and grass) compete for the same food (soil, sunlight and space) and their dynamics are strongly affected by human interactions [Innes_2013, Bauch_2016]. Despite the simplicity of the model, as we show, its nonlinear dynamical behavior is very rich including various types of codimension-one and codimension-two criticalities. We derive analytically the conditions for the onset of these criticalities including the codimension-one transcritical, saddle-node and Andronov-Hopf bifurcations as well as the codimension-two cusp and Bogdanov-Takens bifurcations. We also perform a two-parameter numerical bifurcation analysis to identify far-from-equilibrium criticalities such as limit points of limit cycles that cannot in principle be found analytically. By doing so, we provide a full picture of the mechanisms pertaining to the onset of codimension-one and codimension-two criticalities in a combined way. In the next section, we briefly present the mathematical model. We then provide the analytical results and consequently perform a numerical continuation of the limit cycles and construct the two-parameter bifurcation diagram to fully characterize the different dynamical regimes.
1 Model Description
The model describes the dynamics of a simple mosaic ecosystem of two species, grassland and forest under the human impact. Let and be the proportion of forest and grassland in the ecosystem, with . Humans influence the fate of the evolution by selecting one of the species based on their abundance or lack. Let be the percentage of the human population that prefer forest; thus, is the ratio of people that prefer grassland. Each individual can alter his/her preference to forest or grassland according to some social-depended rules (educational or mimetic), or according to the perceived value of the forest. The coupled model of ordinary differential equations (ODEs) reads [Innes_2013]:
where express the transition rate from forest to grassland due to natural processes and is the rate of transition from grassland to forest mainly induced due to fire modelled by the following relation:
The function is of sigmoid type while the properties of its shape depend from the parameters , and ; gives the maximum value, while , are related to the rate of activation of .
The human influence is represented in the above system of ODEs by the term , where is the potential magnitude of human influence on the ecosystem. The human interaction has a positive or negative feedback on the perceived value of the forest: if the majority prefers forest (i.e. ) the society will enhance forest recreation (i.e ) by reforestation; instead if social opinion favours grassland then the human impact leads to deforestation (i.e ). The perceived value of the forest is given by: .
For more details on the model see Innes et al. 2013[Innes_2013].
2 Analytical Results
We first give a geometrical description for the multiplicity of fixed point solutions. Then, we derive analytically the fixed point solutions and characterize their stability. Consequently, we proceed with the construction of the one-parameter bifurcation diagrams with respect to the human influence as well as the parameter. Then, we provide the conditions for the existence of the saddle-node and Andronov-Hopf bifurcations.
2.1 Fixed Point Solutions-Stability Analysis
Setting the derivative of the second equation of Eq.(1) equal to zero, we get three solutions: , the and .
Starting with and substituting to the first equation of Eq.(1) and setting this to zero, we obtain
The above equation defines a nonlinear function of fixed points for the forest density () with respect to when :
In Fig.1 we have plot for a specific choice of the values of the other parameters of the model.
For , we get a symmetric expression for the parameter for the second branch of fixed points i.e.
This is the mirror case of the fixed points corresponding to with respect to the y-axis (see blue line in Fig.1).
Finally, for and by setting zero the first equation of the system (1), we get:
The above gives rise to another family of solutions of fixed points with given by:
2.1.1 A Geometrical Description of the Multiplicity of Fixed Points
Fig.2 shows the fixed point solutions for , as they result from Eq.(4) and Eq.(5) as intersection points of the graphs of the functions and for various choices of the values of the other parameters. Fig.2(a) depicts the intersection of the graphs for , , , ). Clearly, for there are not fixed points for the branch of solutions. Fig.2(b) depicts the intersection of the graphs for the branch of solutions. As increases, the number of steady state regimes first increases from one to three and then decreases to zero. Finally, it should be considered that the branch of steady state solutions exists for all values of .
For , the intersections of the graphs are depicted in Fig.2(c),(d) as is varied. In particular, in Fig.2(c) we show fixed points corresponding to and in Fig.2(d) the ones corresponding to . As it is shown, for , as is larger than , the number of solutions goes from three to two, indicating that is a bifurcation point. When , there are no solutions corresponding to . For the case of , when , there are three solutions, two of them disappearing at . Again for very high values of , there are no solutions corresponding to .
While this graphical analysis can give a quick idea of the number of solutions for different values of , it cannot give a precise information about the stability and type of bifurcations that may appear. Moreover, dynamic regimes like periodic solutions cannot be detected. Thus, in what follows, we perform a systematic analytical study of the stability of fixed point solutions and their bifurcations.
2.1.2 Stability Analysis of Fixed Points
The stability of the above families of fixed points is determined from the eigenvaluesof the Jacobian matrix of 1 evaluated at the fixed points.
The Jacobian matrix of 1 reads:
Clearly, when or , the Jacobian matrix has an upper triangular form since is vanished. In these cases, the eigenvalues of the Jacobian matrix are real and equal to the diagonal elements. This implies that any loss of hyperbolicity (for or ) cannot be due to Andronov-Hopf bifurcations.
The determinant of the Jacobian matrix for is given by:
Eq.(11) provides the stability relation between the and branches: the stability of one branch implies the instability of the other. Thus, it suffices to study the stability properties of just one of them, for example the solution branch. From Eq.(2), the derivative of with respect to is given by:
Note that both and are always positive (i.e. is a positive increasing function of ). The trace of at is given by:
Thus, for , we examine separately two cases that mark the onset of a critical point: (a) (resulting from (see Eq.(9))) and (b) .
As we prove in the next section, the case (a) corresponds to the appearance of a transcritical bifurcation, while the case (b) corresponds to the appearance of a saddle-node bifurcation.
2.2 Transcritical Bifurcations
Depending on the values of the other parameters, i.e. , , , the value of may change sign.
The upper triangular form of the Jacobian matrix evaluated at implies that , while the second eigenvalue is given by (see Eq.(9)): .
Clearly if is of constant sign around :
then the second eigenvalue changes sign with respect to , when : (when and when ).
If ,i.e. when
then we have a change of stability of solutions
from unstable to stable .
If , i.e. when
the transition is from unstable nodes () (two positives eigenvalues) to saddle nodes (two eigenvalues with opposite sign).
We summarise the above results in the following Theorem.
The point is a critical point for the branch of fixed point solutions defined by . The value of the parameter that defines this criticality is given by . In the case where (16) holds true, there is a change of stability from unstable states when to stable states when . If Eq.(17) holds, then the branch remains locally unstable.
To determine the type of bifurcation at , we proceed as follows.
The Jacobian matrix at reads:
meaning that and the derivative of the system Eq.(1) with respect to the parameter is , meaning that
Then, is a simple bifurcation point [Seydel_2010] and since there are exactly two intersecting branches, we conclude that this type of bifurcation is a transcritical bifurcation.
The stability of the branch can again be determined with the aid of Eq.(7) and Eq.(14). Thus, the solution for the branch can be re-written as:
For , (while if , we have the opposite (non-physical) scenario i.e. ). The determinant along the branch of solutions is:
Note that and assuming that is constant, if , and, if .
Here, we comment that the stability properties of this branch are not independent from the stability of the branch. The trace of the Jacobian at is:
Then, assuming that (and in combination with the sign of the determinant) this is translated as a change from unstable to stable fixed point solutions. But the assumption is identical with Eq.(16), which defines the change of stability on the branch around the critical point , from unstable to stable solutions (see Eq.(16-17)).
To conclude, at the fixed point we have an exchange of stability, which is characteristic of transcritical bifurcation. We state the following theorem:
The system undergoes a transcritical bifurcation at . In the cases where Eq.(16) and hold true, then, as and increase, the two branches and exchange their stability from unstable to stable state and vice versa.
2.3 Saddle-Node and Cusp Bifurcations
2.3.1 Bifurcations on the branch of fixed point solutions
Here, we examine the second case () for the existence of non-hyperbolic points on the branch , (see Eq.(9)).
For this case, we obtain:
The system of Eq.(4),Eq.(24) defines the appearance of critical points (the Jacobian becomes singular). The two nonlinear equations can be solved numerically (e.g. by the Newton-Raphson method) with respect to the unknowns , . We will show analytically, that depending on the value of the parameters , and , the above nonlinear system can have at most two solutions. Furthermore, we show that the corresponding criticalities are saddle-node bifurcations.
From Eq.(24) we see that as , . As then , then , concluding . Fig.3 depicts for different values of the parameter , which confirms the asymptotical behavior of at the boundaries and . In the neighborhood of , it converges to a constant function the , while, in the neighborhood of , it converges asymptotically to . Between the two boundaries, has a bell shape resulting from the first part of , which is . Thus, has two monotonic intervals: the in which it increases and the in which decreases. Both and are positive functions and increases. We study the monotonicity properties of the first part of . We consider the function :
then by setting , Eq.(25) reads
where now . Then, we get the location of the maximum at or equivalently . Substituting this in Eq.(24) we obtain:
Fig.3 shows the (marked with circles) for different values of the parameter . Note that the values of are in a very good agreement with actual maxima of . Depending on the values of , the maximum of can be positive:
Then from the mean value theorem for continuous functions, this implies that the function has two roots correspond to turning points. Since
then in order to have two turning points, a sufficient condition for the parameters , , is given by:
In the critical condition where:
the two turning points collapse to a cusp bifurcation.
In the first case, where the function has two roots, say , the augmented matrix calculated on the roots of reads:
and for , has a rank implying that the the roots of are turning points [Seydel_2010]. The stability properties of the branch near the turning points results in a straightforward manner from the signs of eigenvalues. We assume that the roots of have the following order . This ordering can be achieved if Eq.(16) holds true and simultaneously (which for example is the case with ). Table 1 summarizes the stability along the with respect to the signs of the eigenvalues.
We summarise the above analytical results in the following theorems.
Assuming that the condition (31) holds for the model parameters, then the branch of fixed points posses a cusp bifurcation.
2.3.2 Bifurcations on the branch of fixed point solutions
In this section, we conduct our analysis for the second symmetric branch of fixed point solutions, . Following the previous analysis, we study the transcritical bifurcation at the point and the stability properties of turning points.
The determinant of Jacobian matrix along the is given by:
and similar to the case, the triangular form of the Jacobian matrix gives immediately the eigenvalues on the branch: and . For the simplicity of the presentation, we suppose that hold true the same assumptions for (i.e. Eq.(16)) and for the ordering of the roots of as in the previous case.
The determinant becomes singular at , while if and if . Near the point and as the parameter increases, the branch changes stability from stable to unstable. The critical value of the parameter on the bifurcation point resulting from Eq.(5) by setting is given by: or
The Jacobian matrix evaluated at the point is exactly the same as the one evaluated at (Eq.(18)) and the condition Eq.(19) continues to be valid, concluding that at a transcritical bifurcation appears.
Assuming that , then for . The determinant reads:
concluding that for and for . Since the condition given by Eq.(16) holds true, the sign of the trace of the Jacobian reads:
The above implies that as increases, the stability of the branch changes from unstable to stable .
We summarize the above results with the following theorem.
The stability properties of the branch can easily be calculated from Eq.(11) and from the analysis of the function. Table 2 summarizes the stability along the branch with respect to the sign of the eigenvalues. We observe that for , the branch of solutions remains unstable. At the turning points , the transition is from saddle to unstable nodes (in ) and again back to saddles.
As an example of the above analysis, we present the bifurcation diagram of with respect to for , , (Fig.4). Solid lines depict stable solutions, while dashed lines unstable ones. Note that for this choice of the values of , , , there are no fixed point solutions for (); hence this bifurcation branch corresponds to the locus of steady states with . The horizontal line corresponds to the regime fixed points with . Starting from small values of , the system has two solutions, one stable corresponding to and another unstable for . As increases, the ”S" branch crosses the horizontal line and the two branches exchange their stability through a transcritical bifurcation () at . The ”S" branch also displays two saddle-node bifurcation points () at and () at . Between () and (), the system displays four steady states and just one of these is stable.
|Bifurcation Point||Kind of Bifurcation||Branch|
Returning back to Fig.1, we are now in the position to characterize the types of bifurcations and stability of the solution branches. Fixing the values of the parameters at , , , , the annotated bifurcation diagram of with respect to the parameter is shown in Fig.5.
Solid lines depict stable fixed points, while dashed lines depict unstable ones. As discussed above, there are three branches of steady states: one inverted “S branch corresponding to , an “S" branch corresponding to the locus of steady states with , and the horizontal line corresponding to the fixed points with . In order to understand, from a nonlinear dynamical point of view, the behaviour of the branch of fixed points, we show the bifurcation diagram also for negative values of . In particular, in Fig.5 the stability of each fixed point is reported: are stable nodes with two negative eigenvalues; are saddles with one negative and one positive eigenvalue, and are unstable nodes with two positive eigenvalues. It is apparent, that the branch of the fixed points appears with a mirror “S"-shape with respect to the “‘S" branch (as we proved in the previous section). Both the “S"-shaped branches ( and ) cross the horizontal line giving rise to two transcritical bifurcations ( for the branch, for the branch) on which the branches of fixed points exchange their stability. Thus, looking just at the positive values of , up to the value corresponding to the saddle node bifurcation , there are six steady states: two stable nodes (one for and one for ), three saddles (one for ,one for and one for ) and one unstable node (for ). Multistability and multiplicity of fixed points is observed in a wide range of parameters, whereas as the parameter passes the value of one corresponding to the saddle-node point , only one stable steady state exists.
|Bifurcation Point||Kind of Bifurcation||Branch|
When the human influence is strong enough, the system reaches an balanced equilibrium with 50% of grass and 50% of forest. On the other hand when the human influence is weak, then multiple stable steady states may occur and thus catastrophic transitions are possible between steady states with a very different percent of grass and forest as we discuss in the next section.
2.4 Catastrophic shifts
In an ecosystem, catastrophic shifts occur when the system passes from a stable regime to another completely different regime in an abrupt way as a consequence of environmental or external disturbances. When this happens, the recovery of the previous state is extremely difficult due (in the majority of the cases) to hysteresis effects. From a dynamical point of view, a catastrophic shift may occur in two ways: (a) as a consequence of perturbations of the vector state in a parameter region where the system has multistability, and/or (b) as a consequence of a perturbation of the value of a parameter near criticality.
In the first case, if the system has more than one stable regime for a specific value of the bifurcation parameter (here, the parameter ), then as a consequence of a perturbation in the vector state (here or ), the system may jump from one equilibrium to another. As each stable regime is characterized by its own basin of attraction (the set of all the initial conditions leading to the same stable regime), a perturbation to the vector state may bring it to another basin of attraction which drives the system to a different stable state.
This phenomenon is illustrated in Fig.6, where temporal simulation results are reported for an parameter value belonging to the region of multistability . In this region, there are two stable steady states: one corresponding to and one to . Starting from an initial condition within the basin of attraction of , the system reaches first the steady state and then, when a perturbation is imposed (around t=100), the system finally reaches the steady state (see Fig.6(a),(b)). Indeed, as a consequence of the state vector perturbations, the system trajectory crosses the separatrix of the attraction basins of the two stable states (see Fig.6(c)).
The second way with which a catastrophic transition may occur is when one equilibrium disappears as a consequence of a perturbation of the value of the bifurcation parameter near the critical point. In this case, the system reaches a new stable equilibrium in an abrupt way. Such transitions are observed for example when a parameter crosses a critical value corresponding to a catastrophic bifurcation such as a saddle-node bifurcation point. This case is illustrated in Fig.7, where it is clearly shown that if the system is at the stable steady state just before the critical value of , then a small perturbation of the parameter leads the system to another equilibrium with a sharp transition.
2.5 Andronov-Hopf and Bogdanov-Takens Bifurcations
First, we examine the existence of periodic solutions due to the appearance of Andronov-Hopf bifurcations. Since the two branches , have always real eigenvalues, we seek for Andronov-Hopf points on the third branch which is the . The Jacobian matrix along the branch of solutions reads:
with characteristic polynomial:
If the discriminant of the characteristic polynomial is negative, i.e. and the roots are passing the axis with non zero slope, then the assumptions for the Andronov-Hopf bifurcation are satisfied. In that case, the eigenvalues are summing up to zero or equivalently:
or after short calculations
with roots while must satisfy Eq.(7) and that (taking that ). Eq.(41) is independent from and together with Eq.(7) define a system of equations which has always a solution for each with only restriction the positiveness of the product of . Taking into account Eq.(7) we then obtain:
In the case where , the last inequality holds if . In order to better understand the above behaviour, we performed a two-parameter bifurcation analysis with respect to the parameters , .
2.5.1 Two-parameter Bifurcation Analysis for the Andronov-Hopf Points
where . The function has two monotonic intervals with range ; then if the Eq.(44) has two uneven roots with respect to which can be expressed in an implicit form as:
For the construction of the bifurcation diagrams shown in Fig.8, we used the same values for the parameters , , , resulting to and (these values correspond to the horizontal lines in Fig.8(a)). In the same figure, we also depict the two-parameter continuation of turning points, resulting from the system of equations Eq.(4),(24) and Eq.(5),(24). In Fig.8(b), we depict the continuation of turning points for higher values of and also for . In Fig.8(b), two symmetric branches with respect to the vertical line exist, approaching asymptotically -infinity as (see Fig.8(b)).
An explanation of this asymptotic behavior is provided by using elements of singular two-scale analysis for the function (24). The domain of forest mass variable is . Assuming that , then in the scale of order of i.e. , a bell shape front appears. Since and then one root of (the turning point) belongs to the interval i.e. . As increases, this implies . From Eq. (4) or Eq.(5), we obtain , which is the asymptotic behavior in the neighborhood of .
The second asymptotic behavior results in the slow scale (). Then, the large value of vanishes the first part of concluding . Setting in order to take the second root , we obtain . The corresponding value of is given from Eq.(4) and Eq.(5) resulting to the asymptotic values of i.e.
where the correspond to branches, respectively.
Of particular interest is the case when Eq.(41) holds true along the or branches, respectively (opposite signs in (15)). Based on our analysis, this case occurs when one of the saddle-node points (on branches) collapse with the transcritical points at the branch (which happens when the first root of , ). In that case, the Jacobian matrix has the form:
and immediately we obtain as a double eigenvalue of the matrix. Thus, the critical point for the double zero eigenvalue at the branch can be computed from the system:
For the particular choice of the values of the other parameters, the above system of equations has the solution .
Similarly for , the system is the same, except that and , with .
Setting , the Taylor expansion of the vector field of Eq.(1) around the point reads: