Incommensurate conformable-type three-dimensional Lotka–Volterra model: discretization, stability, and bifurcation

Abstract The classic Lotka–Volterra model is a two-dimensional system of differential equations used to model population dynamics among two-species: a predator and its prey. In this article, we consider a modified three-dimensional fractional-order Lotka–Volterra system that models population dynamics among three-species: a predator, an omnivore and their mutual prey. Biologically speaking, population models with a discrete and continuous structure often provide richer dynamics than either discrete or continuous models, so we first discretize the model while keeping one time-continuous dependent variable in each equation. Then, we analyze the stability and bifurcation near the equilibria. The results demonstrated that the dynamic behaviors of the discretized model are sensitive to the fractional-order parameter and discretization parameter. Finally, numerical simulations are performed to explain and validate the findings, and the maximum Lyapunov exponents is computed to confirm the presence of chaotic behavior in the studied model.


Introduction
The dynamic behaviors studies of nonlinear models emerging in engineering and science play a considerable role in our daily activities. Population dynamics models are widely used in ecology. Understanding the dynamics of prey-predator interactions is one of the most primary operations that shape the framework and function of ecological societies (Rilov, 2009). For instance, some recent experimental studies have demonstrated that ocean acidification influence prey-predator interactions (Kroeker, Sanford, Jellison, & Gaylord, 2014) as ocean acidification makes it difficult for sea snails to burrow from their sea star predators (Jellison, Ninokawa, Hill, Sanford, & Gaylord, 2016). The oldest and most celebrated prey-predator model is the Lotka-Volterra model, independently introduced by Lotka (Lotka, 1925) and Volterra (Volterra, 1926). This model is formulated in ordinary differential equations such that any slight change of the model will lead to a qualitatively different type of behavior. The Lotka-Volterra model may indeed be the simplest possible prey-predator model. Nevertheless, it is a useful tool that contains the essential features of real ecosystems, and serves as a robust foundation from which it is possible to develop more sophisticated generalized models (Korobeinikov & Wake, 1999).
In an ecosystem, we may accept a three-species model that includes prey, predator, and omnivore. An omnivore is a predator that feeds on several trophic levels (Diethelm & Ford, 2002), and it has been seen in food chain models of three or more species. In the present work, we will consider an incommensurate conformable-type fractional-order prey-predator-omnivore model. It is noteworthy that the omnivore considered in the current work is presented as a top scavenger predator, which not only devours the corpses of the predator, but also predates the thoroughbred prey. The investigated system reads as: D a s xðtÞ ¼ rxðtÞð1 À xðtÞÞ À axðtÞyðtÞ À xðtÞzðtÞ, D b s yðtÞ ¼ yðtÞðÀc þ dxðtÞÞ, D c s zðtÞ ¼ xðtÞzðtÞ À bzðtÞ þ eyðtÞzðtÞ À fz 2 ðtÞ, where the conformable fractional derivative Dã s is defined for a function f : ½s, 1Þ ! R, s ! 0, as (Khalil, Al Horani, Yousef, & Sababheh, 2014) and all the constant coefficients a, b, c, d, e, f and r are positive real numbers. In this model, the positive functions xðtÞ, yðtÞ, and z(t) stand for the prey, predator and omnivore densities at time t, respectively. From Equation (2), it has been evinced in (Abdeljawad, 2015) the following necessary fact Dã s f ðtÞ ¼ ðt À sÞ 1Àã f 0 ðtÞ: (3)

Discretization process
Several studies revealed that the discrete-time system exhibits much fruitful dynamic behaviors, such as bifurcations and chaos, than those of its continuous-time system counterpart. Consequently, in this section, we aim to discretize the model (1) making use of the piecewise-constant approximation ) as shown below.
Applying the rule in Equation (3) to the first equation in the system (4), gives the following Bernoulli differential equation: ðt À nhÞ 1Àa dxðtÞ dt þ ðzðnhÞ þ ayðnhÞ À rÞxðtÞ which leads to À x 0 ðtÞ x 2 ðtÞ þ ðr À zðnhÞ À ayðnhÞÞ xðtÞðt À nhÞ 1Àa ¼ r ðt À nhÞ 1Àa : (6) Next, multiply Equation (6) Integrating Equation (7) Let t ! ðn þ 1Þh in Equation (8) and replace x(nh) by x n , yields x nþ1 ¼ x n ðr À z n À ay n Þ rx n þ ðr À z n À ay n À rx n Þ e ÀðrÀz n Àay n Þ h a a : On a similar manner, from the second equation in the system (4), we have dyðtÞ yðtÞ ¼ ðdxðnhÞ À cÞ ðt À nhÞ 1Àb dt: Integrating both sides of Equation (10) on the interval ½nh, tÞ, leads to ln yðtÞ À ln yðnhÞ ¼ ðdxðnhÞ À cÞ For t ! ðn þ 1Þh in Equation (11) and replacing y(nh) with y n , we obtain Again, similar to the first equation, we get z nþ1 ¼ z n ðx n þ ey n À bÞ fz n þ ðx n þ ey n À fz n À bÞ e Àðx n þey n ÀbÞ h c c : As a consequence, the discrete version of the model (4) is given by the following difference equations: x nþ1 ¼ x n ðrÀz n Àay n Þ rx n þðrÀz n Àay n Àrx n Þ e ÀðrÀzn Àayn Þ h a a , y nþ1 ¼ y n e ðdx n ÀcÞ h b b , z nþ1 ¼ z n ðx n þey n ÀbÞ fz n þðx n þey n Àfz n ÀbÞ e Àðxn þeynÀbÞ h c c : 3. Equilibria and analysis of local stability In this section, we explore the local stability for the system of difference Equation (14) around selected equilibria. It is readily verified that the following points are from the equilibria of system (1) Theorem 3.1. The equilibrium E 0 ¼ ð0, 0, 0Þ of system (14) is a saddle point.
Proof. The Jacobian matrix computed at the equilibrium E 0 is given as Thus, E 0 is a saddle point since jk 1 j > 1, jk 2 j < 1 and jk 3 j < 1: Proof. The Jacobian matrix computed at the equilibrium E 1 of the linearization of system (14) is given as and has eigenvalues k 1 ¼ e Àr h a a , which satisfy jk 1 j < 1, k 2 ¼ e ðÀcþdÞ h b b and k 3 ¼ e ð1ÀbÞ h c c : Now, we consider the following cases: (i) if d > c or b < 1, then jk 2 j > 1 or jk 3 j > 1 and E 1 is a saddle point; (ii) if d < c and b > 1, then jk 2 j < 1 and jk 3 j < 1 and E 1 is a sink; (iii) if d ¼ c or b ¼ 1, then jk 2 j ¼ 1 or jk 3 j ¼ 1 and E 1 is a non-hyperbolic.
and r < aðbdÀcÞ eðdÀcÞ : Proof. The Jacobian matrix for system (14) computed at the equilibrium E 2 is given as It's effortless to check that k 2, 3 fulfill the equation Now, jk 1 j < 1 whenever r < aðbdÀcÞ eðdÀcÞ , and according to the Jury conditions (Kot, 2001), jk 2, 3 j < 1 when- (14) is local asymptotically stable if these inequalities are satisfied jC 1 þ C 3 j < 1 þ C 2 , jC 1 À 3C 3 j < 3 À C 2 , and C 2 3 þC 2 ÀC 1 C 3 < 1, where Proof. The Jacobian matrix for system (14) computed at the equilibrium E 3 is given as The characteristic polynomial for the above-mentioned Jacobian matrix is derived as: where the values of C 1 , C 2 , and C 3 are as given in the Equation (16). Now, according to the Jury conditions (Grove & Ladas, 2004), the positive equilibrium E 3 of system (14) is locally asymptotically stable if following conditions are hold: jC 1 þ C 3 j < 1 þ C 2 , jC 1 À 3C 3 j < 3 À C 2 and C 2 3 þ C 2 À C 1 C 3 < 1: w

Neimark-Sacker bifurcation analysis
In this section, we investigate the presence of Neimark-Sacker bifurcation at the positive equilibria E 2 and E 3 of the discrete system (14).
Theorem 4.1. The positive equilibrium point E 2 ¼ ð c d , rðdÀcÞ ad , 0Þ loses its stability via a Neimark-Sacker bifurcation when Hence, the equilibrium E 2 loses its stability via a Neimark-Sacker bifurcation when (14) undergoes a Neimark-Sacker bifurcation for dðb þ rf Þ > cð1 þ rf Þ and aðc À bdÞ > reðc À dÞ if the following conditions are satisfied: where the values of C 1 , C 2 , and C 3 are as given in the Equation (16).
Proof. The characteristic polynomial for the Jacobian matrix of system (14) computed at the equilibrium E 3 is given in Equation (17). According to Lemma 4.1 in (Ali, Saeed, & Din, 2019) for n ¼ 3, the positive equilibrium E 3 of system (14)    The bifurcation diagrams in Figure 2 show that with increasing value of the discretization parameter h with a, b and c values fixed at 0.95 leads to destabilize the equilibrium E 2 through a Neimark-Sacker bifurcation.
Next, we examined the bifurcation and limit cycle for the system (14) at the coexistence state E 3 , choosing r as a bifurcation parameter. Figure 3 demonstrates that E 3 of system (14) is stable for r < 0.275 and loses its stability thought a Neimark-Sacker bifurcation when r ¼ 0.275 and an invariant circle appears for 0:275 < r < 0:595: Thereafter, for r > 0:595 the discrete system (14) return to stabilize. The maximum Lyapunov exponents corresponding to Figures 2(b) and 3(b) are given in Figure 4.

Results and discussion
In this article, the dynamics of the incommensurate conformable-type fractional-order Lotka-Volterra model  with omnivore are investigated. We have shown that the fractional-order parameters a, b and c have impact on the stability of the discrete system. Chaos in the discrete system is studied. It is found that the system exhibits chaotic behavior when increasing b with fixing the fractional-order parameters a and c. The simulation results displayed that the discretization parameter h is a key factor that affects dynamic behavior and is only responsible for a Neimark-Sacker bifurcation and chaos in the dynamic of the system (14). When increasing the discretization parameter h with fixing the fractional-order parameters a, b and c, the discrete conformable-type system is destabilized and chaotic behavior occurs. In the analysis, d and r were chosen as bifurcation parameters for the equilibria E 2 and E 3 , respectively. The results exhibited that the equilibrium E 3 possesses two Neimark-Sacker bifurcation points. To further confirm the chaos, we plot the time series of the system (14), see Figure 5. It is clear that when the bifurcation parameters d and r are increasing, this leads to chaotic behaviour in the system (14).

Disclosure statement
No potential conflict of interest was reported by the authors.