Bifurcation and Stability Analysis of a System of Fractional-Order Di ﬀ erential Equations for a Plant–Herbivore Model with Allee E ﬀ ect

: This article concerns establishing a system of fractional-order di ﬀ erential equations (FDEs) to model a plant–herbivore interaction. Firstly, we show that the model has non-negative solutions, and then we study the existence and stability analysis of the constructed model. To investigate the case according to a low population density of the plant population, we incorporate the Allee function into the model. Considering the center manifold theorem and bifurcation theory, we show that the model shows ﬂip bifurcation. Finally, the simulation results agree with the theoretical studies.


Introduction
Mathematical modeling for various biological problems is considered to be an exciting research area in the discipline of applied mathematics. In the literature, many biological phenomena were modeled and formulated mathematically [1][2][3][4][5][6]. For example, Liu and Xiao established a predator-prey system in discrete time to analyze the local stability and the bifurcation of solutions around the positive equilibrium point [4]. Kangalgil and Kartal analyzed the host-parasite model that led to a system of differential equations of piecewise constant arguments at specific time t, as well as the stability of all obtained equilibrium points; they showed the conditions of flip and Neimark-Sacker bifurcation [6]. Most of these studies are restricted to integer-order differential equations or differential equations with piecewise constant arguments. However, it was seen that many problems in biology, as well as in other fields such as engineering, finance, and economics, could be formulated successfully using fractional-order differential equations [1,5,[7][8][9][10][11][12][13].
The nonlocal property of fractional-order models not only depends on the current state but also depends on its prior historical states [14]. The transformation of an integer-order model into a fractional-order model needs to be precise with respect to the order of differentiation α. However, a small change in α may cause a big change in the behavior of the solutions [15].
Fractional-order differential equations can model complex biological phenomena with non-linear behavior and long-term memory, which cannot be represented mathematically by integer-order differential equations (IDEs) [16,17]. For example, Bozkurt established the glioblastoma multiforme (GBM)-immune system (IS) interaction using a fractional-order differential equation system to include the delay time (memory effect) [5].
In general, the origin of plant-herbivore interactions is derived from predator-prey systems [18,19], which is described in many studies using discrete and continuous time [6,[20][21][22][23]. The usual discrete host-parasite models have the form P n+1 = µP n · f (P n , H n ) H n+1 = cµP n (1 − f (P n , H n )), (1) where P and H are the densities of the host (a plant) and the parasite (a herbivore), µ is the host's inherent rate of increase (µ = e r where r is the intrinsic rate of increase) in the absence of the parasites, c is the biomass conversion constant, and f is the function defining the fractional survival of hosts from parasitism [21]. In Reference [23], the authors considered a system of differential equations of plant-herbivore interactions as follows: K − αx(t)y(t) dy dt = −sy(t) + βx(t)y(t). ( Kartal in Reference [20] generalized the system in Equation (2) by combining discrete and continuous time as follows: where x(t) and y(t) denote the populations of the plant and herbivore, respectively, t is the integer part of t ∈ [0, ∞) and the parameters belonging to R + , r is the growth rate, K is the carrying capacity, and α is the predation rate of the plant species, while s and β represent the death rate and conversion factor of the herbivores, respectively. In our study, we consider a model as a system of FDEs as follows: K − γ f (t)y(t) D α y(t) = β f (t)y(t) − dy(t), (4) where f (t) represents the Holling type II function given by In Equation (4), r is the growth rate of the plant population, K denotes the carrying capacity, γ denotes the predation rate of the plant species, β is the conversion of consumed plant biomass into new herbivore biomass, and d is the per capita rate of death. In Equation (5), e is the encounter rate, which depends on the movement velocity of the herbivore species. The parameter σ (0 < σ ≤ 1) is the fraction of food items encountered that the herbivore ingests, while h is the handling time for each prey item, which incorporates the time required for the digestive tract to handle the item.

Definition 1.
[24] Let f : R + → R be a function, where the fractional integral of order α > 0 is given by provided the right side is pointwise defined on R + .

Definition 2.
[24] Let f : R + → R be a continuous function. The Caputo fractional derivative of order α ∈ (n − 1, n) is given by

Equilibrium Points
Let us consider the system We want to discuss the stability of the system in Equation (8). Let us perturb the equilibrium point by adding ε 1 (t) > 0 and ε 2 (t) > 0, that is Thus, we have and Using the fact f (x, y ) = g(x, y ) = 0, we obtain a linearized system about (x, y) such as where Z = (ε 1 (t), ε 2 (t)), and J is the Jacobian matrix evaluated at the point (x, y), We have B −1 JB = C, where C is a diagonal matrix of J given by where λ 1 and λ 2 are the eigenvalues and B represents the eigenvectors of J. Therefore, we get whose solutions are given by the following Mittag-Leffler functions: and Using the result of Reference [25], if arg(λ 1 ) > απ 2 and arg(λ 2 ) > απ 2 , then η 1 (t) and η 2 (t) are decreasing; consequently, ε 1 (t) and ε 2 (t) are decreasing. Let the solution (ε 1 (t), ε 2 (t)) of Equation (12) exist. If the solution of Equation (12) is increasing, then (x, y) is unstable; otherwise, if (ε 1 (t), ε 2 (t)) is decreasing, then (x, y) is locally asymptotically stable. The equilibrium points of the system in Equation (8) are where β > hd and K > d eσ(β−hd) .

Local Stability
The Jacobian matrix for the system in Equation (8) is given as For Λ 1 = (0, 0), we have the characteristic equation Theorem 1. Let Λ 1 = (0, 0) be the extinction point of the system in Equation (8). Then, the equilibrium point Λ 1 of the system in Equation (8) is a saddle point.
For the case where only the plant population exists, we consider the equilibrium point Λ 2 = (K, 0). The characteristic equation around Λ 2 is as follows: Theorem 2. Assume that Λ 2 = (K, 0) is the equilibrium point of the system in Equation (8). Then, the following statements are true: 1+heσK , then Λ 2 is locally asymptotically stable; (ii) For r > 0, if d < βeσK 1+heσK , then Λ 2 is an unstable saddle point.
then we either attain real or complex conjugates roots with negative real parts, where arg(λ) > απ 2 is equivalent to the Routh-Hurwitz case. Thus, Λ 3 is locally asymptotically stable; then we attain complex conjugates with positive real parts and which implies that Λ 3 is locally asymptotically stable. Proof.
(i) Let us consider the case where ∆ = (A 1 ) 2 − 4A 2 > 0. From where β > . Furthermore, computations show that, for we obtain d > , it is obvious that A 2 > 0. This completes the proof of (i).

Example 1.
In this section, we analyze the stability conditions for the system in Equation (8) that shows a plant-herbivore model with fractional-order differential equations. The values of the parameters are K = 10, eσ = 1.2, h = 1, d = 0.00175, β = 1.5, and γ = 1, where the order of the system is α = 0.9. In Figure 1, we obtain the bifurcation structure of the system in Equation (8) for the initial values (x 0 , y 0 ) = (0.6, 0.75), where the blue graph represents the plant population, and the red graph denotes the herbivore population. Obviously, an increase in the plant population allows the herbivore population to increase. After that, an interaction between both populations occurs. Figure 2 is the per capita for each population, which shows that, after a specific capacity of the habitat, the herbivore population will not be able to find enough food, and this may lead to death or migration. after a specific capacity of the habitat, the herbivore population will not be able to find enough food, and this may lead to death or migration.

Existence and Uniqueness
Considering the system in Equation (8) with initial conditions (0) > 0 and (0) > 0, the initial value problem can be written in the form . Let us assume (0) ≥ and (0) > 0 when > ≥ 0. In this case, the initial value problem can be written as after a specific capacity of the habitat, the herbivore population will not be able to find enough food, and this may lead to death or migration.

Existence and Uniqueness
Considering the system in Equation (8) with initial conditions (0) > 0 and (0) > 0, the initial value problem can be written in the form . Let us assume (0) ≥ and (0) > 0 when > ≥ 0. In this case, the initial value problem can be written as

Existence and Uniqueness
Considering the system in Equation (8) with initial conditions x(0) > 0 and y(0) > 0, the initial value problem can be written in the form . Let us assume x(0) ≥ a and y(0) > 0 when t > σ ≥ 0. In this case, the initial value problem can be written as Definition 3. Assume that C * [0, T] is the class of continuous column vector U(t) whose components is a solution of the initial value problem in Equation (27) if (i) and (ii) hold.
If Definition 4 holds, we obtain the theorem below.

Theorem 4.
Let U ∈ C * [0, T] be a solution of the initial value problem given in Equation (27). Then, U is a unique solution for Equation (27).

Proof. Let us write
Operating with I α , we obtain Now, let F : Then, Therefore, this initial value problem is equivalent to the initial value problem in Equation (27).

Analyzing the Plant-Herbivore Population at Low Density
Allee in 1931 established an important role in the dynamical behavior of populations. He showed that population dynamics with logistic equations in a low population size should be modified with the Allee function in order to represent a realistic phenomenon [26].
By considering the logistic equation, the density increases when the per capita growth rate decreases monotonically; however, it is shown that, in logistic population models with the Allee effect, the per capita growth rate increases to a maximum point at low population density and decreases when the density of the population increases [26]. Many theoretical and laboratory studies showed the essential need of the Allee effect in small populations. Based on the biological studies, the following assumptions are necessary for defining the Allee function: Considering the conditions above, we apply an Allee function at time t to the system in Equation (8) as follows: where t > 0 and (x(0), y(0)) = (x 0 , y 0 ). Moreover, we define E 0 as the Allee coefficient of the plant population, and a(x) = x(t) E 0 +x(t) is the Allee function. The herbivore population is dependent on the plant population. Thus, if then the plant population is not sufficient for the herbivore to exist. For a low population size of the plant population, let us consider the stability conditions around Λ 3 . The Jacobian matrix at Λ 3 is given by Thus, we obtain the characteristic equation where Theorem 5. Let the system in Equation (32) then we either attain real roots or complex conjugates with negative real parts, where arg(λ) > απ 2 is equivalent to the Routh-Hurwitz case, which means that Λ 3 is locally asymptotically stable; then both roots are complex conjugates with positive real parts and which implies that Λ 3 is locally asymptotically stable. Proof.

Theorem 6.
The system in Equation (32) has a unique solution in W = and Proof. Let H ∈ C * [0, T] be a solution of the initial values problem in Equation (32), which holds (t, H(t)) ∈ V, where t ∈ [0, T] and V = [0, T] × W. The Equation (32) can be written as x(t) E 0 +a By operating both sides with I α , we get Let us define the operator F : From Equation (47), we can write which gives Thus, if we choose M such that eσK(β−dh)(rd(1+heσa)−γeσbd)+rd 2 (1+heσa) e 2 σ 2 K(1+heσa)(β−dh) 2 (E 0 +a)M α < 1, we obtain Moreover, we can obtain If we choose N such that . Similarly, one can show that Equation (46) is equivalent to the initial value problem in Equation (32). This completes the proof.

Example 2.
The values of the chosen parameters are similar to those in Example 1. The blue graph represents the plant population, while the red graph denotes the herbivore population. The Allee coefficient is given by E 0 = 0.14. Figure 3 shows the bifurcation structure of the system in Equation (32) under the initial values (x 0 , y 0 ) = (0.6, 0.75). Figure 4 shows the per capita for the plant-herbivore population. We realized that, for a low density of the plant population, the dearth or migration of the herbivore population will occur earlier than expected. The plant population can recover after the herbivore population disappears from that habitat.
Thus, if we choose such that

Example 2.
The values of the chosen parameters are similar to those in Example 1. The blue graph represents the plant population, while the red graph denotes the herbivore population. The Allee coefficient is given by = 0.14. Figure 3 shows the bifurcation structure of the system in Equation (32) Figure 4 shows the per capita for the plant-herbivore population. We realized that, for a low density of the plant population, the dearth or migration of the herbivore population will occur earlier than expected. The plant population can recover after the herbivore population disappears from that habitat.    Moreover, we can obtain If we choose such that Similarly, one can show that Equation (46) is equivalent to the initial value problem in Equation (32). This completes the proof.□

Example 2.
The values of the chosen parameters are similar to those in Example 1. The blue graph represents the plant population, while the red graph denotes the herbivore population. The Allee coefficient is given by = 0.14. Figure 3 shows the bifurcation structure of the system in Equation (32) under the initial values ( , ) = (0.6, 0.75). Figure 4 shows the per capita for the plant-herbivore population. We realized that, for a low density of the plant population, the dearth or migration of the herbivore population will occur earlier than expected. The plant population can recover after the herbivore population disappears from that habitat.

Flip Bifurcation with Discretization Process
In this section, we consider at first the discretization process and the analysis of flip bifurcation. This discretization is an approximation for the right-hand side of the fractional differential equation D α x(t) = f (x(t)), t > 0, where α ∈ (0, 1). We modify our system in Equation (8) in considering the discrete time effect on the model. The discretization of Equation (8) is as follows: The solution of Equation (49) reduces to In repeating the discretization process n times, we get For t ∈ [nh, (n + 1)·h) , while t → (n + 1)·h and α → 1 , we have The Jacobian matrix J of Equation (53) at the equilibrium points is where

Proof. The characteristic equation of J(Λ 3 ) is of the form
where and From Equations (58) and (59), we have Thus, the characteristic equation of J(Λ 3 ) has two eigenvalues, which are If both |λ 1 | < 1 and |λ 2 | < 1, then the equilibrium point is locally asymptotically stable. From we obtain where β < hd(Keσh+1) and Keσh > 1. For h α r 2Γ(α+1) = µ, Equation (61) can be rewritten as follows: where we obtain From we have Considering Equations (63) and (65) together, we obtain This completes the proof.
Thus, we have where According to the flip bifurcation conditions in References [1,32], we obtain that, if δ 5 + δ 3 2 0, the system in Equation (53) undergoes flip bifurcation. Furthermore, if δ 5 + δ 3 2 > 0, then the bifurcation of the second period points is stable, while, for δ 5 + δ 3 2 < 0, it is unstable. According to the flip bifurcation conditions in References [1,32], we obtain tha 0, the system in Equation (53) Figure 6 demonstrate the phase plane portraits of the system in Equation (8) capacity values. Phase plane for = .

Conclusions
In this paper, the biological dynamics of fractional-order differential equatio herbivore model was discussed and analyzed. The local stability of the obtained equ and the existence and uniqueness of the solution in the system in Equation (27) were Example 1). An impressive result was considered in Section 4 where we found that, fo the plant population, the herbivore population disappears. We noticed that the p model is mainly dependent on the plant population size and carrying capacity (see E the other hand, we investigate possible bifurcation types, where we saw that the sys flip bifurcation structure (Section 5). Similar bifurcation studies were observed in herbivore models such as in References [4,17], where they obtained periodic or solutions. Finally, we conclude that the environmental carrying capacity of the plant a strong effect on the system in Equation (27), while the density of the plant species sho effect for the system in Equation (32). The numerical simulations were carried out usin Author Contributions: Bozkurt conceived the study and was in charge of overall directio Bozkurt and Yousef designed the model and set up the main parts of the study. Bozkurt and Y theorems and proved them together. They collected and analyzed the data. Both authors inte and carried out this implementation. Yousef carried the simulations using Matlab 2018. Boz wrote the manuscript and revised it to the submitted form. There was no ghost-writing.
Funding: This research received no external funding. According to the flip bifurcation conditions in References [1,32], we obtain that, if + ≠ 0, the system in Equation (53) undergoes flip bifurcation. Furthermore, if + > 0, then the bifurcation of the second period points is stable, while, for + < 0, it is unstable.   Figure 6 Phase plane for = .

Conclusions
In this paper, the biological dynamics of fractional-order differential equations in a plantherbivore model was discussed and analyzed. The local stability of the obtained equilibrium points and the existence and uniqueness of the solution in the system in Equation (27) were analyzed (see Example 1). An impressive result was considered in Section 4 where we found that, for a low size of the plant population, the herbivore population disappears. We noticed that the plant-herbivore model is mainly dependent on the plant population size and carrying capacity (see Example 2). On the other hand, we investigate possible bifurcation types, where we saw that the system exhibits a flip bifurcation structure (Section 5). Similar bifurcation studies were observed in many plant-

Conclusions
In this paper, the biological dynamics of fractional-order differential equations in a plant-herbivore model was discussed and analyzed. The local stability of the obtained equilibrium points and the existence and uniqueness of the solution in the system in Equation (27) were analyzed (see Example 1). An impressive result was considered in Section 4 where we found that, for a low size of the plant population, the herbivore population disappears. We noticed that the plant-herbivore model is mainly dependent on the plant population size and carrying capacity (see Example 2). On the other hand, we investigate possible bifurcation types, where we saw that the system exhibits a flip bifurcation structure (Section 5). Similar bifurcation studies were observed in many plant-herbivore models such as in References [4,17], where they obtained periodic or quasi-periodic solutions. Finally, we conclude that the environmental carrying capacity of the plant population has a strong effect on the system in Equation (27), while the density of the plant species shows an essential effect for the system in Equation (32). The numerical simulations were carried out using Matlab 2018.
Author Contributions: F.B.Y. conceived the study and was in charge of overall direction and planning. F.B.Y. and A.Y. designed the model and set up the main parts of the study. F.B.Y. and A.Y. set up the theorems and proved them together. They collected and analyzed the data. Both authors interpreted the data and carried out this implementation. A.Y. carried the simulations using Matlab 2018. F.B.Y. and A.Y. wrote the manuscript and revised it to the submitted form. There was no ghost-writing.