Dynamics of a guanaco–sheep competitive system with unilateral and bilateral control

In this paper, based on a guanaco–sheep competitive system, we develop and analyze mathematical models with unilateral and bilateral control for the management of overgrazing. We first analyze the dynamics of the uncontrolled system. It then follows the analysis of the system with impulsive control by differential equation geometry theory. We mainly prove the existence and stability of order-1 periodic solution for unilateral control system and order-2 periodic solution for bilateral control system. Some numerical simulations including the bifurcation diagrams of periodic solution are carried out, which not only verify the validity of the theoretical results, but also reveal some special dynamic phenomena, such as the appearance of higher-order periodic solutions and the existence of parameter intervals with the order change of periodic solution. Comparing the unilateral and bilateral control strategy, we encourage bilateral control rather than unilateral control for the management of sheep species.


Introduction
The guanaco, which distributed in the arid and semiarid areas [1], is a native and wild camelid species in South America [2], and the population size of this species has suffered a sustained decline. Authors of [3] showed that the guanacos occupy only 26% of its original scope. The results of [4] suggested that the range distribution of the guanacos has been reduced by 60% in Argentina, 75% in Chile. This phenomenon is mainly attributed to habitat change and loss, competition with livestock (sheep, cattle and horses) for pasture and the lack of a proper management for livestock [5].
Literature study shows that there is a competitive relationship between the guanaco and sheep species. Puig et al. [6] indicated that competition between guanaco and domestic animals is due to the fact that their diets overlap up to 83%. Baldi et al. [7] and Gordon [8] suggested that guanaco and sheep compete for food and water resources. Ranchers believe that guanaco depletes food supplies, hence reducing productivity of sheep species [9]. So it is very meaningful to investigate the interaction of these two species. In this regard, Laguna et al. [10] developed a three-dimensional mathematical model with two competing herbivores (sheep and guanaco) and one predator, and analyzed the steady state of the occupation of patches for the three species. Daza et al. [11] proposed a delayed differential model capturing the interaction between sheep and guanaco, and investigated coexistence and extinction of these two species.
However, overgrazing or domestic animal overstock has caused vegetation degradation, and the reduction of wildlife is the most serious consequences of vegetation degradation. Pedrana et al. [1] indicated that sheep species seems to have a higher impact on vegetation. Marino et al. [12] showed that grassland degradation and the decline of guanaco species will continue due to overgrazing. Bharucha et al. [13] pointed out that overgrazing modified the plant cover and impaired its productive capacity. So what management programs should be developed by human beings to change this situation and protect guanaco species and grassland? With human intervention, the population size of sheep and guanaco species will change in a short time, and this phenomenon could be well described by impulsive differential equation [14] which has been extensively applied in applications [15][16][17][18][19][20][21][22][23][24][25][26][27]. For example, integrated pest management models were studied in [15][16][17] where the existence of order-1 periodic solution is mainly proved. Meanwhile, the authors of [18] and [19] also considered the optimal control problem. For resource management, Guo et al. [22] analyzed an algae-fish with unilateral control. Moreover, bilateral control strategy is applied in [23,25,27] where the existence and stability of order-2 periodic solution are mainly proved. Huang et al. [27] also analyzed the attraction region of the system, and performed bifurcation diagrams of periodic solution. In addition, bilateral control strategy has also been applied to biodiversity conservation and commodity price adjustment [24,26].
Motivated by the above research, in this paper, by constructing semi-continuous dynamical systems, we develop mathematical models for the management of sheep and guanaco species. On the one hand, the number of sheep species will be limited to avoid overgrazing; on the other hand, the economic interests of farmers should also be protected while keeping sheep from overgrazing. In practical operation, the number of sheep is more controllable compared with guanaco, and its number monitoring is more convenient and feasible, so we take the number of sheep as the state variable of feedback control in the design of state feedback control strategy. Based on these two aspects, we consider unilateral and bilateral control strategies to keep the sheep species in a suitable range.
The rest of the paper is organized as follows. In Sect. 2, we develop guanaco-sheep competitive sys-tems with unilateral and bilateral control strategy. In Sect. 3, the existence of the order-1 and order-2 periodic orbit is proved by successor function. And stability of the order-1 and order-2 periodic orbit is verified by analogue of the Poincaré criterion. In Sect. 4, the numerical simulations are carried out with a special case to validate the theoretical results. Finally, a brief conclusion is presented in Sect. 5.

Model formulation
For the guanaco and sheep species, their competitive ability is related to their body size. Tilman et al. [28] considered that the competitiveness of these two competing herbivores is not in the same hierarchy. Daza et al. [11] also pointed out that there exists hierarchical competition between guanaco and sheep, the guanaco is assumed to be the superior species, and sheep is the inferior one. So we assume that the effect of inferior species on superior species satisfies Holling II functional response cy 1+y [29][30][31], and the influence of higher competitors on lower competitors satisfies a linear function ex. Based on the classical two-species competitive system [32], we get a modified competitive system as follows: where x(t), y(t), respectively, stands for the density of guanaco and sheep species. a, d, respectively, represents the intrinsic increasing rate of the guanaco and sheep species; c measures the strength of the effect of sheep species on the guanaco; e measures the strength of the effect of the guanaco on the sheep species; b, f , respectively, stands for the intra-specific competition coefficient of the guanaco and sheep species.
As the application of semi-continuous dynamic system in overgrazing management, we mainly investigate guanaco-sheep competitive systems with unilateral and bilateral control strategies. In practical operation, the number of sheep is more controllable compared with guanaco, and its number monitoring is more convenient and feasible, so we take the number of sheep as the state variable of feedback control in the design of state feedback control strategy. On the one hand, desertification caused by overgrazing should be gradually mitigated; on the other hand, the interests of farmers should also be protected. So we are committed to keeping the number of sheep within a reasonable range [h 1 , h 2 ]. To this end, we extend system (1) by incorporating two control thresholds y = h 1 and y = h 2 . When the number of sheep species decreases to h 1 , herdsman can increase the number of sheep by a certain proportion, while reducing the number of guanacos by a certain amount. When the quantity of sheep species reaches to h 2 , people can release guanaco and harvest sheep suitably. Thus, two guanaco-sheep competitive systems with unilateral control are established as follows: and When two kinds of control are considered simultaneously, we have the following competitive system with bilateral control: All the parameters in systems (2), (3) and (4) are positive, and 0 < α 1 , α 2 < 1.

Main results
In this section, we mainly discuss the dynamical properties of system (1) and the existence and stability of periodic solution of systems (2), (3) and (4).

Dynamics of free system (1)
For system (1), it is easy to determine that trivial equilibrium O(0, 0) is always unstable. The equilibrium otherwise it is unstable. Positive equilibrium (x * , y * ) satisfies the following equation: After simple calculation, we have x * = d− f y * e , and y * satisfies the following second-order algebraic equation: Obviously, the curve l 2 := {(x, y)|w(x, y) = 0} has a horizontal asymptote y = −1 and a vertical asymptote y = a−c b , respectively. Combining with differential equation geometry theory, we analyze the existence and stability of the positive equilibrium of system (1) from the following two cases. Case 1. If a > c, then vertical asymptote y = a−c b lies on the right of y-axis. The location of the two isoclines l 1 and l 2 is shown in Fig. 1. Fig. 1a). The Jacobian matrix of system (1) evaluated at E * 1 is So E * 1 is locally stable. By constructing Dulac function D(x, y) = 1 xy , it is easy to find that system (1) has no limit cycle. So we conclude that E * 1 is globally stable.
is a saddle point.

Case 2.
If a < c, then vertical asymptote y = a−c b lies on the left of y-axis. The location of the two isoclines l 1 and l 2 is shown in Fig. 2.  Fig. 2b). It is easy to determine that E * is a saddle point.

Existence and stability of the order-1 periodic solution
In this work, we mainly prove the existence and stability of periodic solution when the unique positive equilibrium E * 1 (x * 1 , y * 1 ) is unstable. It then follows that the proof the existence of order-1 periodic solution of system (2).

Existence of the order-1 periodic solution of system (2)
As is shown in Fig. 3, the dashed lines indicate the separatrix of saddle point which divide the first quadrant into four regions. Line l 3 intersects with the phase set . We now illustrate the existence of the order-1 periodic solution of system (2).
Proof Denote the intersection point of the curve l 2 and the line M 1 by B 1 . According to the phase diagram of system (2), there must exist a point A 1 (x A 1 , h 1 ) ∈ N 1 such that the trajectory from A 1 intersects with M 1 at point B 1 , then jumps to A + 1 due to the impulsive effect.
Thus, we only need to find another point A 2 ∈ N 1 such that G(A 2 ) > 0. Consider the intersection point A 2 of the dashed line l 3 and N 1 , the trajectory from A 2 intersects with N 1 at point B 2 , and then jumps to A + 2 due to the impulsive effect. Based on the condition Fig. 4. According to the continuity of successor functions, there must exist a point A ∈ N 1 between A 1 and A 2 satisfying that G(A) = 0, i.e., the trajectory AB ∪ B A forms an order-1 periodic solution of system (2).
It then follows the proof of the stability of order-1 periodic solution of system (2).

Stability of order-1 periodic solution of system
(2)
Remark 1 According to Theorems 1 and 2, we can see that when the lower threshold of the number of sheep raised by herdsmen is limited and the state feedback control is carried out, through the timely rounding up of guanaco (τ 1 ) and the supply of the number of sheep (α 1 ) raised, we can ensure that the number of the two species presents a stable periodic change and avoid the situation that the number of one species is too small and even extinct.

Existence and stability of the order-1 periodic solution of system (3)
Next, we mainly investigate the existence of the order-1 periodic solution of system (3).

Fig. 4
The existence of order-1 periodic solution of system (3) As is shown in Fig. 4, line l 4 intersects with the phase set Proof Select a point A 3 (x A 3 , (1 + α 2 )h 2 ) in phase set N 2 where x A 3 is sufficiently small, the trajectory from A 3 intersects with M 2 at point B 3 , and then jumps to A + 3 due to the impulsive effect. Since x A 3 is sufficiently small, then the successor function of point A 3 satisfies that G( Next, we only need to find another point A 4 ∈ N 2 such that G(A 4 ) < 0. Consider the intersection point A 4 of the dashed line l 4 and N 2 , the trajectory from A 4 intersects with M 2 at point B 4 , and then jumps to A + 4 due to the impulsive effect. Since Fig. 4. By the continuity of successor function, there must exist a point A between A 3 and A 4 satisfying that G(A) = 0, i.e., there exists an order-1 periodic solution of system (3).

Remark 2
The proof of the stability of order-1 periodic solution of system (3) is similar to the proof of theorem 2, so we omit it here.
Remark 3 According to Theorem 1 and Remark 2, we can see that when the upper threshold of the number of sheep is limited and the state feedback control is carried out, through the timely supply of guanaco (τ 2 ) and reasonable harvest of sheep (α 2 ), we can also ensure that the number of the two species presents a stable periodic change.

Existence and stability of the order-2 periodic solution
In this subsection, we mainly prove the existence and stability of order-2 periodic solution of system (4).

Existence of the order-2 periodic solution of system (4)
Next, we will prove the existence of order-2 periodic solution of system (4). For convenience, we suppose that the isocline l 1 intersects with impulsive set M 2 and phase set N 1 at points G 1 (x G 1 , h 2 ) and G 2 (x G 2 , (1 + α 1 )h 1 ), respectively. Isocline l 2 and the line l 3 intersect with impulsive set M 1 at point h 1 ), respectively. l 5 and l 6 intersect with N 1 and N 2 at point , respectively. Besides, there must exist a point S(x S , (1 − α 2 )h 2 ) ∈ N 2 , such that the trajectory from point S intersects with impulsive set M 1 at point D 3 .

then system (4) has an order-2 periodic solution.
Proof Choose a point A 1 ∈ N 1 , where A 1 is very closed to y-axis, i.e., x A 1 is positive and small enough. The trajectory starting from point A 1 will intersect with the impulsive set M 2 at point B 1 (x B 1 , h 2 ). Based on the second impulsive function, point B 1 is mapped to point The trajectory starting from point C 1 will intersect with the impulsive set M 1 at point D 1 (x D 1 , h 1 ), then D 1 is on the right side of the line l 3 . Based on the first impulsive function, point D 1 is mapped to point A + 1 ∈ N 1 , where point A + 1 is the order-2 successor point of A 1 , and In addition, x A 1 is positive and small enough, so we have Fig. 5a). Choose another point A 2 ∈ N 1 , and A 2 is sufficiently closed to point F 2 and it is on the left side of point F 2 . Similar to the discussion above, the trajectory starting from point A 2 intersects with M 2 at point By the continuity of successor function, there is a point A between A 1 and A 2 satisfying G(A) = 0, then AB ∪ BC ∪ C D ∪ D A forms an order-2 periodic solution, as is shown in Fig. 5b.
It then follows that the stability of order-2 periodic solution of system (4).

Stability of the order-2 periodic solution
Let T (t) = (ξ(t), η(t)) be the order-2 periodic solution of system (4). Define Fig. 5b). Next, we prove the stability of the order-2 periodic solution.
Remark 4 According to Theorems 4 and 5, we can see that for the bilateral control that both the upper and lower thresholds of the number of sheep are limited and the state feedback control is carried out, through reasonable adjustment of the harvest and supply strategy of the two populations (α 1 , α 2 , τ 1 , τ 2 ), their number will maintain a sustained and stable periodic change.

Examples and numerical simulations
In this section, some numerical simulations will be presented to verify the theoretical results in the previous section. Different types of periodic solutions are performed with different control parameters. Example 1 Select model parameters as follows: a = 0.45, b = 0.7, c = 0.9, d = 0.4, e = 1, f = 0.3. If we just take unilateral control measure, then we have following two guanaco-sheep competitive sys-tems with unilateral control: and In order to get the order-1 periodic solution of system (6), we select control parameters h 1 = 0.15, α 1 = 0.25 and τ 1 = 0.02. The trajectory starting from the initial value (x 0 , y 0 ) = (0.4, 0.5) will tend to be an order-1 periodic solution (see Fig. 6). Figure 6d is a partial enlargement of Fig. 6c.
Keep other parameters unchanged but τ 1 = 0.2, then the trajectory from the same initial value also tends to be an order-1 periodic solution (see Fig. 7). Obviously, comparing Figs. 6d and 7c , we find that system (6) exhibits different types of order-1 solution with different values of τ 1 .
Keep other parameters unchanged, choosing the control parameter τ 2 = 0.05, and the trajectory from same initial value also tends to be an order-1 periodic solution (see Fig. 9). It is obvious that system (7) exhibits different types of order-1 solutions with different values of τ 2 by comparing Figs. 8c and 9c.
Example 2 In order to illustrate the existence of order-2 periodic solution of system (4) in Sect. 3.3, by using the same model parameters, we get the following specific system To get the order-2 periodic solution of system (8), we select control parameters h 1 = 0.25, h 2 = 0.75, α 1 = 3) tends to be an order-2 periodic solution (see Fig. 10). Keep other parameters unchanged but τ 1 = 0.075, and it is easy to find that the trajectory from the same initial value tends to be an order-2 periodic solution (see Fig. 11).
Example 3 To show the bifurcation diagram of system (4), we choose τ 2 ∈ (0, 0.9) as bifurcation parameter and h 1 = 0.25, h 2 = 0.75, α 1 = 0.6, α 2 = 0.45, τ 1 = 0.075. From Fig. 12, we can see that system (4) can exhibit order 1, 2, 3 and 4 periodic orbits with τ 2 ∈ (0, 0.9). More precisely, there exist an order-1 periodic orbit for τ 2 ∈ (0, 0.258), and then an order-2 periodic solution bifurcates from the order-1 periodic solution at τ 2 ≈ 0.258. Also, system (4) exists order-2 periodic orbit for τ 2 ∈ (0.258, 0.45) and order-3 periodic solution for τ 2 ∈ (0.45, 0.792). Similarly, system (4) appears order-3 periodic solution at τ 2 ≈ 0.45 and order-4 periodic solution at τ 2 ≈ 0.792. Choose τ 2 = 0.79, 0.8 and keep other parameters fixed, we show the existence of order-3 and order-4 periodic solution of system (4) (see Fig. 13). Obviously, system (4) has different types of periodic solutions with the change of parameter τ 2 . And if we select τ 1 , α 1 or α 2 as bifur- Remark 5 Numerical results show that there is always a periodic solution when the conditions of the theorems in Sect. 3 are satisfied. In addition, we also note that for the two models of unilateral control (4) and (3), the first pays more attention to the minimum economic interests of herdsmen, and emphasizes that the number of sheep should not be lower than an economic threshold in state control. The second focuses on the maximum benefits that herdsmen can obtain. At this time, the number of guanacos in the environment is at a relatively low level. Corresponding to them, after adopting the bilateral control model (4), on the one hand, the interests of herdsmen have been fully guaranteed; on the other hand, the number of guanacos in the environment has also been maintained at a relatively reasonable level.

Conclusions
In this paper, we develop and analyze a mathematical model for the management of overgrazing based on a guanaco-sheep competitive system with unilateral and bilateral control. The case of that the unique positive equilibrium E * 1 of the basic ODE model is unstable is mainly discussed. More precisely, the existence of order-1 and order-2 periodic solution is verified by the method of successor function. Also, the stability of order-1 and order-2 periodic solution is proved by applying the geometry theory and analogue of Poincaré criterion.
Numerical results show that if the control parameters α 1 , τ 1 are given suitably, then system (2) exists an order-1 periodic solution. It means that guanaco and sheep species can be maintained periodically in certain level. If we change the value of τ 1 and keep other parameters unchanged, then system (2) also has an order-1 periodic solution. Besides, system (2) exhibits different types of order-1 periodic solution with different values of τ 1 . There exist similar results for system (3). In addition, if the control parameters α 1 , τ 1 , α 2 , τ 2 are chosen suitably, then system (4) has an order-2 periodic solution. In the research of mathematical model, the exploration of the influence of different parameters on the system has always been highly valued by scholars [33][34][35][36][37][38]. By drawing bifurcation diagram of periodic solution with τ 2 ∈ (0, 0.9), it is easy to find that system (4) exists high-order periodic solution, such as order-3 and order-4 periodic solutions. Compared with the results obtained in [24,32], the proposed control models in this paper exhibit different dynamical behaviors, such as the different types of order-1 and order-2 periodic solutions, high-order periodic solutions and the existence of parameter intervals with the order change of periodic solution.
The research shows that if unilateral control measure is adopted, then system (2) or (3) has order-1 periodic solution, it means that population size of these two species changes periodically. For system (2), the aim of this control measure is to avoid the extinction of the sheep species and guarantee the farmers' profits, but it may not alleviate desertification caused by overgraz- ing. For system (3), this unilateral control strategy is mainly to reduce the competitive pressure of the guanaco from the sheep species and may alleviate desertification caused by overgrazing; it will be beneficial to prevent the reduction of wildlife. Combining these two aspects, then bilateral control strategy is adopted in system (4), this control strategy can not only effectively prevent the loss of farmers' interests due to the low number of sheep species, but also prevent the guanaco from being subjected to great competitive pressure due to the large number of sheep species, and the desertification caused by overgrazing may be mitigated. Comparing the unilateral control and bilateral control strategies, we encourage bilateral control rather than unilateral control for the management of sheep species.
sive mapping φ 2 maps it to the phase point Q + 2 . Q + 2 is called the order-1 successor point of P 2 , and G(P 1 ) = y Q + 1 − y P 1 is called the successor function of point P 2 . If Q + 2 coincides with P 2 , then the trajectory P 2 Q 2 ∪ Q 2 P 2 forms an order-1 periodic solution, as is shown in Fig. 14a. 3. For a certain point P 1 ∈ N 2 , the trajectory from P 2 intersects with M 1 at point Q 1 , then the impulsive mapping φ 1 maps it to the phase point Q + 1 ∈ N 1 . And, the trajectory from Q + 1 intersects with M 2 at point S 1 , then the impulsive mapping φ 2 maps it to the phase point S + 1 ∈ N 2 . S + 1 is called order-2 successor point of P 1 . If S + 1 coincides with P 1 , then the trajectory P 1 Q 1 ∪ Q 1 Q + 1 ∪ Q + 1 S 1 ∪ S 1 P 1 forms an order-2 periodic solution, as is shown in Fig. 14b. Besides, point S + 1 is called the order-2 successor point of P 1 , and G(P 1 ) = y S + 1 − y P 1 is called the order-2 successor function of point P 2 .