# Insight into Caputo Fractional-Order Extension of Lotka–Volterra Model with Emphasis on Immigration Effect

Authors: A. El-Mesady; Omar Bazighifan; Serkan Araci
Journal: Journal of Mathematics
 (2023)

---

## Abstract

For over a hundred years, the model of prey–predator has been handled extensively. It has been proved that Lotka–Volterra (LV) models are not stable from all previous studies. In these models, there is divergent extinction of one species or cyclic oscillation. There is no solution for asymptotic stability for the prey–predator models. In this paper, a Caputo fractional-order one prey-two predators’ model with an immigration effect is investigated. Then, the effect of adding a small immigration term to the prey or predators’ populations is studied. The extension is both the inclusion of a second predator and the use of Caputo derivatives. Although LV systems have been extensively investigated, there are no previous studies on our proposed nonlinear modifications. The Grünwald–Letnikov method was applied to find numerical solutions for the proposed model to validate our findings. Several different cases were studied to reach confirmed conclusions. The effect of the variation of fractional-order derivative is handled in this study. From all the results, the effect of adding the small immigration terms is positive for stability. Hence, it can be concluded that the proposed small immigration terms can stabilize the natural populations of the one prey-two predators’ model. An optimal control strategy for the proposed model is shown at the end of the paper.


---

## Body

## 1. Introduction

Ecology is the scientific study that relies on the distribution and abundance of living things and how these factors are affected by interactions between the environment and these living things, see [1]. Ecology relies heavily on several sciences, especially pedology, geography, chemistry, geology, physics, and meteorology, see for example [2]. According to Xia and Chen [3], geology, insolation, climate, and the living things that have the same habitat are the essential physical properties of the environment of the living things. Ecology may also be considered as a biological study that can also include studies in different fields such as cells and proteins in cellular biology, nucleic acids in biochemistry, and molecular biology. Ecology can include other disciplines such as zoology, botany, population science, ecosystems, and community studies [4]. Many biomathematical models, depending on laboratory experiments and statistical data [5], have been proposed to describe nerve systems, viruses, and deoxyribonucleic acid (DNA). By finding the analytical and numerical solutions of such biomathematical models, we can get a clear understanding of the underlying phenomena. Then, the control strategies can be designed.The interaction between the prey and predators is the main feature that reveals the behavior of the prey and predators. There is a mathematical model formulated by Lotka [6] and Volterra [7] that describes the interaction between the prey and predator. Dynamics of predator–prey systems in arthropods were studied in [8]. In [9], the authors presented the beginnings and development of the predator–prey theory. A prey–predator model with a prey refuge and a Holling type III response function was studied in [10]. The dynamical behavior of an exploited system consisting of two preys and a predator, that is, being harvested, is investigated in [11]. The impact of immigration and refuge on the dynamics of a three-species system in which one predator feeds on one of two competing species is examined in [12]. The authors in [13] examined the predator–prey-pathogen relationship, in which the predator affects the infection’s transmission rate in its prey. A predator–prey-pathogen model in which the predator is a natural specialist and infected prey can enter refugia of a fixed size to avoid predator attack was proposed in [14]. The bifurcation analysis of a predator–prey model with prey refuge is discussed in [15]. Recently, the classical Lotka–Volterra (LV) model has been modified by several researchers by studying several factors such as disease, Allee effects, and fear mechanisms. For the classical LV model, the mathematical properties indicate either divergent extinction or cyclic oscillation of one species, see [16]. The extinction of the prey in any closed LV system will eventually lead to the death of the predators. Hence, the persistent prey–predator systems will exhibit cyclic oscillations if there are no additional stabilizing mechanisms. However, in wild prey–predator systems, stable coexistence has been observed in [17]. Nevertheless, the lynx-hare interaction case [18] seems to be a unique natural prey–predator system. The coexistence of spider wasps (the predator) and spiders (the prey) in [17] is considered an example of the stable coexistence of both prey and predator populations. Generally, in natural prey–predator systems, an additional stabilizing mechanism should be available.In recent times, fractional calculus is considered a very important topic since it has an important memory effect. Several fields, such as biology, finance, economics, engineering, and medicine, have used fractional calculus to solve their problems; see for example [19] and the references therein. There exist several definitions for fractional-order derivatives. The most popular definition of fractional-order derivative is the Caputo-type derivative [20]. Recently, many research papers have been dedicated to fractional-order systems [21]. Owolabi and Shikongo [22] introduced the fractional operator method for a multimutation and intrinsic resistance model. A mathematical model for multimutation and drug resistance with fractional derivative was introduced by [23]. In [24], the authors handled a fractional-order model and studied the stability using the Routh-Hurwitz criteria. Sene and Srivastava in [25] presented a new stability notion of the fractional differential equations with exogenous input and proposed the generalized Mittag-Leffler input stability conditions. Rashid et al. [26] proved a multiple fixed-point result for partially ordered s-distance spaces under θ,ϕ,ψ-type weak contractive condition. In [27], a fractional-order prey–predator model with discontinuous harvest and functional response of the Holling III type is considered. Kumar et al. [28] provided a comparative study for a fractional LV model in the Caputo sense using Haar wavelet and Adams-Bashforth-Moulton methods. The homotopy perturbation method (HPM), an analytical method for nonlinear problems, was used in [29] to show how the Lotka-Volterra equations of fractional-order temporal derivatives can be solved. With the use of the Lotka–Volterra prey–predator model and logistic growth of prey species, a discrete fractional-order model was introduced in [30]. Several stability characteristics of the states, including Mittag-Leffler stability, practical stability, and stability with regard to sets, were addressed via impulsive control techniques in [31]. For a good survey on the LV model, see [32–35]. Tahara et al. [36] studied the effect of immigration on a Lotka-Volterra model. However, no work has been performed on the Caputo fractional-order Lotka–Volterra model with the immigration effect. Our proposed model is an extension to the model in [36], where the proposed model is a one prey-two predators’ model.The aim of this work and the contributions that have been made are as follows. We develop a Caputo fractional-order one prey-two predators’ model with an immigration effect in this paper. Then, we study the effect of adding a small immigration term to the prey or predator populations. The extension is both the inclusion of a second predator and the use of Caputo derivatives. The reason for selecting this particular Caputo definition is that, unlike the Riemann–Liouville fractional derivative, the Caputo fractional derivative allows for the inclusion of initial conditions in the model. The Caputo fractional derivative of a constant is also zero, in contrast to the Riemann–Liouville fractional derivative. The proposed model, along with the initial conditions, will be the main model that is covered in detail in this study. Although LV systems have been extensively investigated, there are no previous studies on our proposed nonlinear modifications. This study comes to the conclusion that very little immigration into either the prey population or the predators’ population stabilizes the Lotka-Volterra system. All fluctuations in both the prey and predator populations will be averaged out by adding positive immigration. Also take note of the fact that an increase in immigration alone can alter the population’s qualitative dynamics. The following few points highlight the proposed research’s novelty:(1)
A newly developed fractional model, defined in the Caputo sense, for examining the new dynamics of the one prey-two predators’ model with the effect of immigration.(2)
The stability analysis is presented.(3)
To validate the theoretical facets and results of the suggested model, an effective numerical technique is adopted.(4)
The obtained results demonstrate the effectiveness of the suggested model in providing some new insights into the dynamical behavior of the proposed model.(5)
An optimal control strategy is discussed.This paper is divided as follows: Section2 presents the fundamentals and basic definitions that will be utilized later in other sections. Section 3 displays the fractional-order Lotka-Volterra model with immigration effect. Non-negativity, existence, and uniqueness are discussed in Section 4. Stability analysis is introduced in Section 5. We introduced the numerical simulations based on the Grünwald–Letnikov method in Section 6. An optimal control strategy is presented in Section 7. Finally, Section 8 has the conclusion of the paper.


## 2. Fundamentals and Definitions

Some of the fundamental definitions that will be utilized later in other sections are provided in this section, see [37, 38].Definition 1.
The Riemann–Liouville fractional integral for a continuous functionf:0,+∞⟶ℝ is given by the following equation:(1)Itα0RLft=1Γα∫0tt−xα−1fxdx,α∈0,1,t>0.Definition 2.
The Riemann–Liouville fractional derivative for a continuous functionf:0,+∞⟶ℝ can be represented by the following equation:(2)Dtα0RLft=1Γ1−αddt∫0tt−x−αfxdx,α∈0,1,t>0.Definition 3.
The Caputo fractional derivative for a continuous functionf:0,+∞⟶ℝ is defined by the following equation:(3)Dtα0Cft=1Γ1−α∫0tt−x−αddxfxdx,α∈0,1,t>0.
Grünwald–Letnikov (GL) definition of fractional derivatives is given by the following equation:(4)Dtαaft=limh⟶0h−α∑i=0t−a/h−1iαift−ih,where [.] refers to the integer part.
By applying relation (5) formulated from relation (4), we can calculate the fractional-order derivatives numerically. The numerical approximation of α-th derivative at the points kh, (k = 1, 2, ...) can be written as follows [39]:(5)Dtkαk−Lm/hft≈h−α∑i=0k−1iαiftk−i=h−α∑i=0kciαftk−i,where Lm refers to the length of memory, tk=kh,h is step size, and ciα (i = 0, 1, …, k are binomial coefficients. The following expression can be used to calculate ciα as follows:(6)c0α=1,ciα=1−1+αici−1α.
The nonlinear fractional differential equation in the formDtαaft=vxt,t has a general numerical solution that can be expressed by the relations (5) and (6) as follows:(7)ftk=vftk,tkhα−∑i=pkciαftk−i.
A short memory principle can be used. Then, the lower index of the sums in relations (7) will be p = 1 for k<Lm/h and p=k−Lm/h for k>Lm/h. Scherer et al. [40] studied the existence of bounded unique solutions for the Grünwald–Letnikov definition of fractional derivatives. A complete study of the stability of the Grünwald–Letnikov method has been introduced in [41]. Petrás [42] proposed the numerical schemes that we have applied to our model to get the graphics.


## 3. Proposed Fractional Model

Here, we propose the Caputo fractional Lotka–Volterra model with few prey and predators immigrants to reflect the effect of the memory. Also, we study the effect of adding an immigration termDx to the prey population or adding immigration terms GyandJz to the predators populations.The Caputo fractional one prey-two predators model with few immigrants can be represented as follows:(8)Dtα0Cx=ax−bxy−cxz+Dx,Dtα0Cy=exy−fy+Gy,Dtα0Cz=hxz−iz+Jz,where a,b,c,e,f,h, and i are positive real constants. The parameter a refers to the reproduction rate of the prey, b refers to the rate at which the first predator consumes the prey, and c refers to the rate at which the second predator consumes the prey. The parameter e refers to the birth rate of the first predator for each prey captured, and f refers to the mortality rate of the first predator. The parameter h refers to the birth rate of the second predator for each prey captured, and i refers to the mortality rate of the second predator. The prey population is represented by x, the first predator population is represented by y, and the second predator population is represented by z.Dtα0C refers to the Caputo fractional derivative of order α. The proposed model is a generalization to the model introduced in [36]. The immigration function can be represented in three ways as follows:(9)Dx=d,d≥0,dx,d>0,(10)Gy=g,g≥0,gy,g>0,(11)Jz=j,j≥0,jz,j>0,where d represents the number of prey immigrants, g represents the number of first predator immigrants, and j represents the number of second predator immigrants. Also, d/x represents the proportion of prey immigrants, g/y represents the proportion of first predator immigrants, and j/z represents the proportion of second predator immigrants. We proposed the following assumptions for equations (9)–(11): x≠0 when Dx=d/x, y≠0 when Gy=g/y, and z≠0 when Jz=j/z.In this paper, to study the effect of the immigration factor on the Caputo fractional one prey-two predators model (system (8)), we consider the following six cases:Case A1: immigrants of the prey (Dx=d,Gy=0,Jz=0)Case B1: immigrants of the first predator (Dx=0,Gy=g,Jz=0)Case C1: immigrants of the second predator (Dx=0,Gy=0,Jz=j)Case D1: few immigrants of the prey (Dx=d/x,Gy=0,Jz=0)Case E1: few immigrants of the first predator (Dx=0,Gy=g/y,Jz=0)Case F1: few immigrants of the second predator (Dx=0,Gy=0,Jz=j/z)We also handle the following six cases of the immigration factor:Case A2: migrants of the prey (i.e.,Dx=−d,Gy=0,Jz=0)Case B2: migrants of the first predator (i.e.,Dx=0,Gy=−g,Jz=0)Case C2: migrants of the second predator (i.e.,Dx=0,Gy=0,Jz=−j)Case D2: few migrants of the prey (i.e.,Dx=−d/x,Gy=0,Jz=0)Case E2: few migrants of the first predator (i.e.,Dx=0,Gy=−g/y,Jz=0)Case F2: few migrants of the second predator (i.e.,Dx=0,Gy=0,Jz=−j/z)


## 4. Existence, Uniqueness, and Non-Negativity

We established several findings regarding the existence and uniqueness of solutions for system (8) using fixed point theory findings, such as the Schauder and Banach fixed point theorems. We follow [43] in this section. The following compensation is used for the right-hand sides of system (8):(12)ψ1t,x,y,z=ax−bxy−cxz+Dx,ψ2t,x,y,z=exy−fy+Gy,ψ3t,x,y,z=hxz−iz+Jz.This yields the following system:(13)Dtα0Cx=ψ1t,x,y,z,Dtα0Cy=ψ2t,x,y,z,Dtα0Cz=ψ3t,x,y,z,x0=x0,y0=y0,z0=z0.Clearly,ψ1,ψ2,ψ3∈C0,T×ℝ3⟶ℝ upon using integral operator ItαC on both sides of equation (13) and using the initial conditions, one has the following equation:(14)xt=x0+∫0tt−ωα−1Γαψ1ω,xω,yω,zωdω,yt=y0+∫0tt−ωα−1Γαψ2ω,xω,yω,zωdω,zt=z0+∫0tt−ωα−1Γαψ3ω,xω,yω,zωdω.LetB=C0,τ×C0,τ×C0,τ be Banach space equipped with the norm x,y,z=supt∈0,τx+y+z. We defined an operator Ξ=Ξ1,Ξ2,Ξ3:B⟶B by the following equation:(15)Ξ=Ξ1x,y,zΞ2x,y,zΞ3x,y,zt=x0+∫0tt−ωα−1Γαψ1ω,xω,yω,zωdωy0+∫0tt−ωα−1Γαψ2ω,xω,yω,zωdωz0+∫0tt−ωα−1Γαψ3ω,xω,yω,zωdω.We consider the following assumptions to establish the results:There are real positive valuesAψ1,Bψ1,Cψ1,Dψ1,Aψ2,Bψ2,Cψ2,Dψ2,Aψ3,Bψ3,Cψ3,Dψ3>0, such that(16)ψ1t,x,y,z≤Aψ1x+Bψ1y+Cψ1z+Dψ1,ψ2t,x,y,z≤Aψ2x+Bψ2y+Cψ2z+Dψ2,ψ3t,x,y,z≤Aψ3x+Bψ3y+Cψ3z+Dψ3.(ii)
There are constantsEψ1,Eψ2,Eψ3>0, such that for x,y,z,x¯,y¯,z¯∈B, we have the following equation:(17)ψ1t,x,y,z−ψ1t,x¯,y¯,z¯≤Eψ1x−x¯+y−y¯+z−z¯,ψ2t,x,y,z−ψ2t,x¯,y¯,z¯≤Eψ2x−x¯+y−y¯+z−z¯,ψ3t,x,y,z−ψ3t,x¯,y¯,z¯≤Eψ3x−x¯+y−y¯+z−z¯.Theorem 1.
The considered system (8) has at least a solution under assumption (i).Proof.
SupposeV=x,y,z∈B:x,y,z≤k⊂B is a convex and closed set, then to prove that Ξ:V⟶V is continuous, suppose ∀x,y,z∈B, we have the following equation:(18)Ξ1x,y,z=supt∈0,τx0+∫0tt−ωα−1Γαψ1ω,xω,yω,zωdω≤x0+1Γα∫0tt−ωα−1ψ1ω,xω,yω,zωdω≤ταAψ1+Bψ1+Cψ1k+Dψ1Γ1+α,and similarily(19)Ξ2x,y,z≤ταAψ2+Bψ2+Cψ2k+Dψ2Γ1+α,Ξ3x,y,z≤ταAψ3+Bψ3+Cψ3k+Dψ3Γ1+α.
This implies thatΞV⊂V. Also ψ1,ψ2,ψ3 are continuous so is Ξ. For t1<t2∈0,τ, we have the following equation:(20)Ξ1x,y,zt2−Ξ1x,y,zt1=∫0t2t2−ωα−1Γαψ1ω,xω,yω,zωdω−∫0t1t1−ωα−1Γαψ1ω,xω,yω,zωdω≤Aψ1+Bψ1+Cψ1k+Dψ1Γ1+αt2α−t1α.
It is clear that the right side of equation (20) goes to zero as t2⟶t1 and similarly we obtain the following equation:(21)Ξ2x,y,zt2−Ξ2x,y,zt1=∫0t2t2−ωα−1Γαψ2ω,xω,yω,zωdω−∫0t1t1−ωα−1Γαψ2ω,xω,yω,zωdω≤Aψ2+Bψ2+Cψ2k+Dψ2Γ1+αt2α−t1α,Ξ3x,y,zt2−Ξ3x,y,zt1=∫0t2t2−ωα−1Γαψ3ω,xω,yω,zωdω−∫0t1t1−ωα−1Γαψ3ω,xω,yω,zωdω≤Aψ3+Bψ3+Cψ3k+Dψ3Γ1+αt2α−t1α.
Consequently,
Ξ1x,y,zt2−Ξ1x,y,zt1+Ξ2x,y,zt2−Ξ2x,y,zt1+Ξ3x,y,zt2−Ξ3x,y,zt1⟶0 as t2⟶t1. Hence, we have Ξx,y,zt2−Ξx,y,zt1⟶0 as t2⟶t1.
Therefore, we can conclude thatΞ is equi-continues and bounded. Hence, Ξ is uniformly continuous and completely continuous based on the theorem of Arzelá-Ascoli. From the theorem of Schauder’s on fixed point, the operator Ξ has at least a fixed point. Hence, system (8) has at least a solution.Theorem 2.
IfΩ=maxEψ1τα/Γ1+α,Eψ2τα/Γ1+α,Eψ3τα/Γ1+α<1 and the assumption (ii) is verified, then system (8) has a unique solution.Proof.
Supposex,y,z,x¯,y¯,z¯∈B, then consider(22)Ξ1t,x,y,z−Ξ1t,x¯,y¯,z¯≤Eψ1ταΓ1+αx−x¯+y−y¯+z−z¯,Ξ2t,x,y,z−Ξ2t,x¯,y¯,z¯≤Eψ2ταΓ1+αx−x¯+y−y¯+z−z¯,Ξ3t,x,y,z−Ξ3t,x¯,y¯,z¯≤Eψ3ταΓ1+αx−x¯+y−y¯+z−z¯.
Then, from equation (22), we obtain Ξx,y,z−Ξx¯,y¯,z¯≤Ωx−x¯+y−y¯+z−z¯. Hence Ξ is a contraction operator and the system under study (8) has a unique solution.Lemma 1 (See [44]).
Let0<α≤1. Suppose that vt∈Ca,b and Dtαt0Cvt∈Ca,b. If Dtαt0Cvt≥0,t∈a,b, then vt is a nondecreasing function for each t∈a,b. If Dtαt0Cvt≤0,t∈a,b, then vt is a nonincreasing function for each t∈a,b.Only non-negative solutions are of importance to us because of the biological significance. The following theorem proves that the solutions of model (8) are non-negative. Suppose ℝ+ is the set of all non-negative real numbers and Λ+=x,y,z∈Λ:x,y,z∈ℝ+.Theorem 3.
For the model (8), all the solutions that start in Λ+ are non-negative.Proof.
Firstly, we prove that the solutionsxt which start in Λ+ are non-negative, that is, xt≥0 for all t≥t0. Suppose that is not true, then there exists a constant t1≥t0 such that(23)xt>0,t0≤t<t1,xt1=0,xt1+<0.
Based on equation (23) and the first equation of system (8), we have the following equation:(24)Dt1αt0Cxt1|xt1=0≥0.
According to Lemma1, we have xt1+≥0, which contradicts the fact xt1+<0. Therefore, we have xt≥0 for all t≥t0. Using the same manner as above, we have yt≥0,zt≥0 for all t≥t0.


## 5. Stability Analysis

We can find the equilibrium points for the proposed model (8) by solving the following system:(25)ax−bxy−cxz+Dx=0,exy−fy+Gy=0,hxz−iz+Jz=0.The formulation of the previous system (25) is based on the values of the functions Dx and Gy, and Tables 1–4 show the equilibrium points and the corresponding eigen values for all the proposed cases of Dx, Gy, and Jz. Now, we discuss the stability of the equilibrium points as follows:(i)
If the eigenvalue is purely imaginary, the system will oscillate for all time.(ii)
A complex eigenvalue with negative real parts means you have an oscillation of decreasing amplitude until eventually you reach zero.(iii)
Also, for the complex eigen values, when the real part is positive, the system is unstable and behaves as an unstable oscillator.(iv)
A purely negative real eigenvalue means that the solutions are exponential and decay directly to zero.(v)
When the eigenvalue is real, positive, the system is unstable.In Tables1–4, EigJ refers to the eigen values of the Jacobian matrix J that can easily be computed by solving J−λI=0 for λ. We have used the symbol EigJ with some cases since the corresponding eigen values in these cases have a complicated form.Lemma 2 (see [45]).
Consider the following fractional-order system:(26)Dtαt0Cx=fx,x0=x0,0<α≤1,x∈ℝn.As stated above, the equilibrium points of system (8) can be found by solving the following equation: fx=0. These points are locally asymptotically stable if all eigenvalues λi of the Jacobian matrix J=∂f/∂x evaluated at the equilibrium points satisfy the Matignon conditions:(27)argλi>απ2.The stable and unstable regions for0<α<1 are shown in Figure 1. It is clear that there are many eigen values in Tables 1–4. Hence, in what follows and based on Lemma 2, we will discuss the stability of some equilibrium points in Tables 1–4.(i)
From Table1, the corresponding eigen values for the trivial equilibrium point 0,0,0 are λ1=a,λ2=−f,andλ3=−i. Hence, argλ1=0<απ/2 for all 0<α<1,argλ2=π>απ/2 for all 0<α<1, and argλ3=π>απ/2 for all 0<α<1. Consequently, the trivial equilibrium point 0,0,0 is unstable.(ii)
From Table1, the corresponding eigen values for the equilibrium point f/e,a/b,0 are λ1=fh−ei/e,λ2=−1af,andλ3=−−1af. If fh>ei, then argλ1=0<απ/2 for all 0<α<1, argλ2=π/2, and argλ3=π/2. Consequently, the equilibrium point f/e,a/b,0 is unstable.(iii)
From Table1, the corresponding eigen values for the equilibrium point 0,0,j/i are λ1=−f,λ2=−i,andλ3=ai−cj/i. If ai<cj, then argλ1=π<απ/2 for all 0<α<1,argλ2=π, and argλ3=π<απ/2 for all 0<α<1. Consequently, the equilibrium point 0,0,j/i is asymptotically stable.(iv)
From Table4, the corresponding eigen values for the equilibrium point 0,0,−−1j/i are λ1=−f,λ2=−2i,andλ3=a+c−1j/i. Hence, argλ1=π>απ/2 for all 0<α<1,argλ2=π>απ/2, and argλ3=arctancj/ai<απ/2for all. Consequently, the equilibrium point 0,0,−−1j/i is unstable.The stability of the remaining equilibrium points can be discussed similarly. Also, Lyapunov stability test can be applied. Lyapunov functions for systems of fractional differential equations in biology have been constructed in literature, see [46] and the references therein.Figure 1
Stability region of the fractional-order system (8), where 0<α<1.Table 1
Equilibrium points and eigen values for the original case, case A1, case B1, and case C1.
CaseEquilibrium pointsEigen valuesOriginal case0,0,0a,−f,−if/e,a/b,0±−1af,fh−ei/ei/h,0,a/c±−1ai,−f+ei/hCase A1−d/a,0,0a,−de/a−f,−dh/a−if/e,de+af/bf,0−de/2f±b2e2f2d2e2−4def2−4af3/2bef2fh/e−ii/h,0,dh+ai/ci−f+ei/h−dh/2i±c2h2i2d2h2−4dhi2−4ai3/2chi2Case B1af−bg/ae,a/b,0−b2eg±b2e2−4a3f+4a2bg+b2g2/2abe, afh−bgh−aei/ae0,g/f,0−f,a−bg/f,−ii/h,gh/fh−ei,a+bgh/−fh+ei/cEig0−bi/h−ci/hegh/fh−ei−f+ei/h0ha+bgh/−fh+ei/c00Case C1ai−cj/ah,0,a/c−afh−aei+cej/ah−c2h2j±−4a3c2h4i+4a2c3h4j+c4h4j2/2ach20,0,j/i−f,−i,ai−cj/if/e,a+cej/fh−ei/b,ej/−fh+eiEig0−bf/e−cf/eea+cej/fh−ei/b00ehj/−fh+ei0fh/e−iTable 2
Equilibrium points and eigen values for case D1, case E1, and case F1.
CaseEquilibrium pointsEigen valuesCase D1−1d/a,0,02a,−−1de/a−f,−−1dh/a−i−−1d/a,0,02a,−1de/a−f,−1dh/a−if/e,a+de2/f2/b,0−bde3f2±b2e2f4d2e4−de2f3−af5/bef4fh/e−ii/h,0,a+dh2/i2/c−f+ei/h−cdh3i2±c2h2i4d2h4−dh2i3−ai5/chi4Case E1f−b2g/a2/e,a/b,0−b3eg±b2e2−a5f+a3b2g+b4g2/a2befh−b2gh/a2/e−i0,−g/f,0−2f,a+bg/f,−i0,g/f,0−2f,a−bg/f,−ii/h,−gh/fh−ei,a+bgh/fh−ei/cEig0−bi/h−ci/h−egh/fh−ei−2f+2ei/h0ha+bgh/fh−ei/c00i/h,gh/fh−ei,a−bgh/fh−ei/cEig0−bi/h−ci/hegh/fh−ei−2f+2ei/h0ha−bgh/fh−ei/c00Case F1i−c2j/a2/h,0,a/c−f+ei−c2j/a2/h−c3h2j±c2h4−a5i+a3c2j+c4j2/a2ch20,0,−j/i−f,−2i,a+cj/i0,0,j/i−f,−2i,a−cj/if/e,a+cej/−fh+ei/b,−ej/−fh+eiEig0−bf/e−cf/eea+cej/−fh+ei/b00−ehj/−fh+ei02fh/e−2if/e,a−cej/−fh+ei/b,ej/−fh+eiEig0−bf/e−cf/eea−cej/−fh+ei/b00ehj/−fh+ei02fh/e−2iTable 3
Equilibrium points and eigen values for case A2, case B2, and case C2.
CaseEquilibrium pointsEigen valuesCase A2d/a,0,0a,de/a+f,dh/a−if/e,−de+af/bf,0de/2f±b2e2f2d2e2+4def2−4af3/2bef2, fh/e−ii/h,0,−dh+ai/ci−f+ei/h, dh/2i±c2h2i2d2h2+4dhi2−4ai3/2chi2Case B2af+bg/ae,a/b,0b2eg±b2e2−4a3f−4a2bg+b2g2/2abeafh+bgh−aei/ae0,−g/f,0−f,a+bg/f−ii/h,gh/−fh+ei,a+bgh/fh−ei/cEig0−bi/h−ci/hegh/−fh+ei−f+ei/h0ha+bgh/fh−ei/c00Case C2ai+cj/ah,0,a/c−afh−aei−cej/ahc2h2j±−4a3c2h4i−4a2c3h4j+c4h4j2/2ach20,0,−j/i−f,−i,ai+cj/if/e,a+cej/−fh+ei/b,ej/fh−eiEig0−bf/e−cf/eea+cej/−fh+ei/b00ehj/fh−ei0fh/e−iTable 4
Equilibrium points and eigen values for case D2, case E2, and case F2.
CaseEquilibrium pointsEigen valuesCase D2−d/a,0,02a,−de/a−f,−dh/a−id/a,0,02a,de/a−f,dh/a−if/e,a−de2/f2/b,0bde3f2±b2e2f4d2e4+de2f3−af5/bef4, fh/e−ii/h,0,a−dh2/i2/c−f+ei/h, cdh3i2±c2h2i4d2h4+dh2i3−ai5/chi4Case E2f+b2g/a2/e,a/b,0b3eg±b2e2−a5f−a3b2g+b4g2/a2be, fh+b2gh/a2/e−i0,−−1g/f,0−2f,a+−1bg/f,−i0,−1g/f,0−2f,a−−1bg/f,−ii/h,−gh/−fh+ei,a+bgh/−fh+ei/cEig0−bi/h−ci/h−egh/−fh+ei−2f+2ei/h0ha+bgh/−fh+ei/c00i/h,gh/−fh+ei,a−bgh/−fh+ei/cEig0−bi/h−ci/hegh/−fh+ei−2f+2ei/h0ha−bgh/−fh+ei/c00Case F2i+c2j/a2/h,0,a/c−f+ei+c2j/a2/h, c3h2j±c2h4−a5i−a3c2j+c4j2/a2ch20,0,−−1j/i−f,−2i,a+c−1j/i0,0,−1j/i−f,−2i,a−c−1j/if/e,a+cej/fh−ei/b,−ej/fh−eiEig0−bf/e−cf/eea+cej/fh−ei/b00−ehj/fh−ei02fh/e−2if/e,a−cej/fh−ei/b,ej/fh−eiEig0−bf/e−cf/eea−cej/fh−ei/b00ehj/fh−ei02fh/e−2i


## 6. Numerical Simulation

### 6.1. Analysis of Results

Through this paper, we substitutex0=5,y0=5,andz0=10. A time simulation of the proposed model is presented by applying the Grünwald–Letnikov method. First, we show the effect of changing the fractional order on the one prey-two predator model without immigration (original case) (i.e., Dx=0,Gy=0,Jz=0.) By taking a=b=c=0.1,e=h=0.3, and f=i=0.2, Figure 2 depicts how the population dynamics of the one prey-two predators system are affected by changing the value of the fractional-order α. The phase space of the one prey-two predators system is shown in Figure 3.Figure 2
Population dynamics of the one prey-two predators system withDx=0,Gy=0,andJz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1; and 3D-plot of xt,yt, and zt for α=1.Figure 3
Phase space for the one prey-two predators system, whereDx=0,Gy=0,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1.For Case A1, we haveDx=d,Gy=0,andJz=0. By taking a=b=c=0.1,d=0.01,e=h=0.3, and f=i=0.2, Figure 4 depicts how the population dynamics are affected by changing the value of the fractional-order α. The phase space of the one prey-two predators is shown in Figure 5.Figure 4
Population dynamics of the one prey-two predators system withDx=d=0.01,Gy=0,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, α=0.75,0.85,0.9,1, and 3D-plot of xt,yt, and zt for α=1.Figure 5
Phase space for the one prey-two predators system, whereDx=d=0.01,Gy=0,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1.For Case C1, we haveDx=d,Gy=0,andJz=j. By taking a=b=c=0.1,j=0.01,e=h=0.3, and f=i=0.2, Figure 6 depicts how the population dynamics of the one prey-two predators system are affected by changing the value of the fractional-order α. The phase space of the one prey-two predators system is shown in Figure 7.Figure 6
Population dynamics of the one prey-two predators system withDx=0,Gy=0,Jz=j=0.01,a=b=c=0.1,e=h=0.3,f=i=0.2,α=0.75,0.85,0.9,1, and 3D-plot of xt,yt, and zt for α=1.Figure 7
Phase space for the one prey-two predators system, whereDx=0,Gy=0,Jz=j=0.01,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1.For the Case E1, we haveDx=0,Gy=g/y,andJz=0. By taking a=b=c=0.1,g=0.01,e=h=0.3, and f=i=0.2, Figure 8 depicts how the population dynamics of the one prey-two predators system are affected by changing the value of the fractional-order α. The phase space of the prey-two predators system is shown in Figure 9.Figure 8
Population dynamics of the one prey-two predators system withDx=0,Gy=g/y=0.01/y,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2,α=0.75,0.85,0.9,1, and 3D-plot of xt,yt, and zt for α=1.Figure 9
Phase space for the one prey-two predators system, whereDx=0,Gy=g/y=0.01/y,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1.For the Case B2, we haveDx=0,Gy=−g,andJz=0. By taking a=b=c=0.1,g=0.01,e=h=0.3, and f=i=0.2, Figure 10 depicts how the population dynamics of the one prey-two predators system are affected by changing the value of the fractional-order α. The phase space of the one prey-two predators system is shown in Figure 11.Figure 10
Population dynamics of the one prey-two predators system withDx=0,Gy=−g=−0.01,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2,α=0.75,0.85,0.9,1, and 3D-plot of xt,yt, and zt for α=1.Figure 11
Phase space for the one prey-two predators system, whereDx=0,Gy=−g=−0.01,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, α=0.75,0.85,0.9,1.For the Case F2, we haveDx=0,Gy=0,andJz=−j/z=−0.01/z. By taking a=b=c=0.1,j=0.01,e=h=0.3, and f=i=0.2, Figure 12 depicts how the population dynamics of the one prey-two predators system are affected by changing the value of the fractional-order α. The phase space of the one prey-two predators system is shown in Figure 13.Figure 12
Population dynamics of the one prey-two predators system withDx=0,Gy=0,Jz=−j/z=−0.01/z,a=b=c=0.1,e=h=0.3,f=i=0.2,α=0.75,0.85,0.9,1, and 3D-plot of xt,yt, and zt for α=1.Figure 13
Phase space for the one prey-two predators system, whereDx=0,Gy=0,Jz=−j/z=−0.01/z,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1.


### 6.2. Discussion of Results

The periodic orbital relationships between the populations of predator and prey cannot be abolished over time, as demonstrated by the classical LV systems. Keep in mind that if the prey goes extinct, the predator population will also go extinct. This does not, however, imply that if the predator population goes extinct, the prey will follow suit. Instead, the population of prey will continue to exist over time and even increase. There is no alternative stable convergence in the classical models of LV, and the only other solution is the cyclic oscillation for the LV models. However, predator–prey relationships in the wild frequently exhibit steady coexistence of both predator and prey populations. The current research demonstrates that extremely few migrants, whether they are predators or prey, produce stability in the LV systems.The biological meaning of immigration changes whend (or g or j) and d/x (or g/y or j/z) are introduced. A relatively small number of immigrants are continuously added to the population, which may occur if the environment appeals to newcomers based on habitat quality. A very tiny amount of immigration into the populations of prey or predators stabilizes the LV systems, see Figures 2–13. All fluctuations in the population of predators and prey will be averaged out by the addition of positive immigration. Additionally, the findings show that a positive immigration component is sufficient to alter the predator–prey model’s population dynamics. According to this study, cyclic populations can be stabilized by introducing a little amount of immigration. Over time, there are usually at least a few immigrants in most natural communities. These few migrant species are enough to keep the prey–predator models stable. This could explain some of the stability in natural prey–predator models that has been seen. Also, Figures 2–13 demonstrate how the system is dependent on the system’s history and how sensitive it is to changes in the fractional-order α. The observed behavior from the numerical simulations shows that the states have a decreasing behavior as time is increased while the fractional-order α is decreased. The effect of the fractional derivative on the behavior of the dynamics of the proposed model is noticed to show its stabilization effect.


## 6.1. Analysis of Results

Through this paper, we substitutex0=5,y0=5,andz0=10. A time simulation of the proposed model is presented by applying the Grünwald–Letnikov method. First, we show the effect of changing the fractional order on the one prey-two predator model without immigration (original case) (i.e., Dx=0,Gy=0,Jz=0.) By taking a=b=c=0.1,e=h=0.3, and f=i=0.2, Figure 2 depicts how the population dynamics of the one prey-two predators system are affected by changing the value of the fractional-order α. The phase space of the one prey-two predators system is shown in Figure 3.Figure 2
Population dynamics of the one prey-two predators system withDx=0,Gy=0,andJz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1; and 3D-plot of xt,yt, and zt for α=1.Figure 3
Phase space for the one prey-two predators system, whereDx=0,Gy=0,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1.For Case A1, we haveDx=d,Gy=0,andJz=0. By taking a=b=c=0.1,d=0.01,e=h=0.3, and f=i=0.2, Figure 4 depicts how the population dynamics are affected by changing the value of the fractional-order α. The phase space of the one prey-two predators is shown in Figure 5.Figure 4
Population dynamics of the one prey-two predators system withDx=d=0.01,Gy=0,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, α=0.75,0.85,0.9,1, and 3D-plot of xt,yt, and zt for α=1.Figure 5
Phase space for the one prey-two predators system, whereDx=d=0.01,Gy=0,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1.For Case C1, we haveDx=d,Gy=0,andJz=j. By taking a=b=c=0.1,j=0.01,e=h=0.3, and f=i=0.2, Figure 6 depicts how the population dynamics of the one prey-two predators system are affected by changing the value of the fractional-order α. The phase space of the one prey-two predators system is shown in Figure 7.Figure 6
Population dynamics of the one prey-two predators system withDx=0,Gy=0,Jz=j=0.01,a=b=c=0.1,e=h=0.3,f=i=0.2,α=0.75,0.85,0.9,1, and 3D-plot of xt,yt, and zt for α=1.Figure 7
Phase space for the one prey-two predators system, whereDx=0,Gy=0,Jz=j=0.01,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1.For the Case E1, we haveDx=0,Gy=g/y,andJz=0. By taking a=b=c=0.1,g=0.01,e=h=0.3, and f=i=0.2, Figure 8 depicts how the population dynamics of the one prey-two predators system are affected by changing the value of the fractional-order α. The phase space of the prey-two predators system is shown in Figure 9.Figure 8
Population dynamics of the one prey-two predators system withDx=0,Gy=g/y=0.01/y,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2,α=0.75,0.85,0.9,1, and 3D-plot of xt,yt, and zt for α=1.Figure 9
Phase space for the one prey-two predators system, whereDx=0,Gy=g/y=0.01/y,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1.For the Case B2, we haveDx=0,Gy=−g,andJz=0. By taking a=b=c=0.1,g=0.01,e=h=0.3, and f=i=0.2, Figure 10 depicts how the population dynamics of the one prey-two predators system are affected by changing the value of the fractional-order α. The phase space of the one prey-two predators system is shown in Figure 11.Figure 10
Population dynamics of the one prey-two predators system withDx=0,Gy=−g=−0.01,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2,α=0.75,0.85,0.9,1, and 3D-plot of xt,yt, and zt for α=1.Figure 11
Phase space for the one prey-two predators system, whereDx=0,Gy=−g=−0.01,Jz=0,a=b=c=0.1,e=h=0.3,f=i=0.2, α=0.75,0.85,0.9,1.For the Case F2, we haveDx=0,Gy=0,andJz=−j/z=−0.01/z. By taking a=b=c=0.1,j=0.01,e=h=0.3, and f=i=0.2, Figure 12 depicts how the population dynamics of the one prey-two predators system are affected by changing the value of the fractional-order α. The phase space of the one prey-two predators system is shown in Figure 13.Figure 12
Population dynamics of the one prey-two predators system withDx=0,Gy=0,Jz=−j/z=−0.01/z,a=b=c=0.1,e=h=0.3,f=i=0.2,α=0.75,0.85,0.9,1, and 3D-plot of xt,yt, and zt for α=1.Figure 13
Phase space for the one prey-two predators system, whereDx=0,Gy=0,Jz=−j/z=−0.01/z,a=b=c=0.1,e=h=0.3,f=i=0.2, and α=0.75,0.85,0.9,1.


## 6.2. Discussion of Results

The periodic orbital relationships between the populations of predator and prey cannot be abolished over time, as demonstrated by the classical LV systems. Keep in mind that if the prey goes extinct, the predator population will also go extinct. This does not, however, imply that if the predator population goes extinct, the prey will follow suit. Instead, the population of prey will continue to exist over time and even increase. There is no alternative stable convergence in the classical models of LV, and the only other solution is the cyclic oscillation for the LV models. However, predator–prey relationships in the wild frequently exhibit steady coexistence of both predator and prey populations. The current research demonstrates that extremely few migrants, whether they are predators or prey, produce stability in the LV systems.The biological meaning of immigration changes whend (or g or j) and d/x (or g/y or j/z) are introduced. A relatively small number of immigrants are continuously added to the population, which may occur if the environment appeals to newcomers based on habitat quality. A very tiny amount of immigration into the populations of prey or predators stabilizes the LV systems, see Figures 2–13. All fluctuations in the population of predators and prey will be averaged out by the addition of positive immigration. Additionally, the findings show that a positive immigration component is sufficient to alter the predator–prey model’s population dynamics. According to this study, cyclic populations can be stabilized by introducing a little amount of immigration. Over time, there are usually at least a few immigrants in most natural communities. These few migrant species are enough to keep the prey–predator models stable. This could explain some of the stability in natural prey–predator models that has been seen. Also, Figures 2–13 demonstrate how the system is dependent on the system’s history and how sensitive it is to changes in the fractional-order α. The observed behavior from the numerical simulations shows that the states have a decreasing behavior as time is increased while the fractional-order α is decreased. The effect of the fractional derivative on the behavior of the dynamics of the proposed model is noticed to show its stabilization effect.


## 7. Optimal Control Strategy

Here, we proposed three time-dependent control techniques in the same time. We calculated the optimal values foru1,u2,andu3 that would reduce the oscillations in x,y,z. Additionally, the costs are minimized. For more illustrations for different optimal control strategies, see [47–49] and the references therein. The proposed system for the optimal control is as (28)Dtα0Cx=ax−bxy−cxz+u1x,Dtα0Cy=exy−fy−u2y,Dtα0Cz=hxz−iz−u3z.ForU=u1t,u2t,u3t∈U, where(29)U=U∈L∞0,T3|0≤ui≤uf,0≤uf≤1,i=1,2,3,the objective functional is defined by the following equation:(30)JU∗=minU∈UJU,where(31)JU=∫0TA1yt+A2zt−A3xt+1/2∑i=13Ci⋅ui2tdt.Here,A1 represents the weight for the number of the first predator, A2 represents the weight for the number of the second predator, A3 represents the weight for the number of the prey, and Cii=1,2,3 stand for weight coefficients and express the unit costs of the controls at different levels. We will determine the required conditions for the optimal control U∗ based on Pontryagin’s maximum principle [50].Theorem 4.
Letx,y,z be the state solutions of problem (27) and (28) and u1∗,u2∗,u3∗ be the optimal controls, then there are adjoint functions satisfying the following equation:(32)Dtα0Cλ1=∂ϕ∂x=−A3+a−by−cz+u1λ1+eyλ2+hzλ3,Dtα0Cλ2=∂ϕ∂y=A1+−bxλ1+ex−f−u2λ2,Dtα0Cλ3=∂ϕ∂z=A2+−cxλ1+hx−i−u3λ3,with the transversality conditions λiT=0andi=1,2,3. Furthermore, the optimal controls u1∗,u2∗,u3∗ are as follows:(33)u1∗=min1,max0,−λ1xC1,u2∗=min1,max0,λ2yC2,u3∗=min1,max0,λ3zC3.Proof.
Using Pontryagin’s maximum principle, ifU∗∈U is considered an optimal solution (27) and (28), then there exists a continuous mapping λ:0,T⟶ℝ, where λt=λ1t,λ2t,λ3t is the adjoint vector, which is obtained from the following equation:(34)Dtα0Cx=∂ϕ∂λ1,Dtα0Cy=∂ϕ∂λ2,Dtα0Cz=∂ϕ∂λ3,that is,(35)Dtα0Cλ1=∂ϕ∂x,Dtα0Cλ2=∂ϕ∂y,Dtα0Cλ3=∂ϕ∂z.
The Hamiltonian functionϕ is as follows:(36)ϕ=A1yt+A2zt−A3xt+12∑i=13Ci.ui2t+λ1.Dtα0Cx+λ1.Dtα0Cy+λ3.Dtα0Cz.
The adjoint system can be obtained as a result. To obtain the characterization of the optimal controlU∗, we use the following equation:(37)∂ϕ∂u1=∂ϕ∂u2=∂ϕ∂u3=0.
Hence, we obtain the optimal control characteristic (29).
In the remaining part of this section, we depict how the population dynamics of the prey–predators system are affected on the basis of controlling functions. Figure14 shows the effect of applying the control signals u1,u2,andu3 on the state xt. The solid lines show the state xt without applying the control signals u1,u2,andu3. The dashed lines show the state xt after applying the control signals u1,u2,andu3. Figure 15 shows the effect of applying the control signals u1,u2,andu3 on the state yt. The solid lines show the state yt without applying the control signals u1,u2,andu3. The dashed lines show the state zt after applying the control signals u1,u2,andu3. Figure 16 shows the effect of applying the control signals u1,u2,andu3 on the state zt. The solid lines show the state zt without applying the control signals u1,u2,andu3. The dashed lines show the state zt after applying the control signals u1,u2,andu3. It is clear that applying the control signals reduces the oscillations. Figure 17 shows the control profile u1,u2,andu3 at different values of α=1,0.9, and 0.8.Figure 14
Dynamics ofxt where a=0.1,b=0.1,c=0.1,e=0.3,f=0.2,h=0.3,i=0.2,C1=0.0002,C2=0.0005,C3=0.0003,A1=0.01,A2=0.02,andA3=0.03.Figure 15
Dynamics ofzt where a=0.1,b=0.1,c=0.1,e=0.3,f=0.2,h=0.3,i=0.2,C1=0.0002,C2=0.0005,C3=0.0003,A1=0.01,A2=0.02,andA3=0.03.Figure 16
Dynamics ofzt where a=0.1,b=0.1,c=0.1,e=0.3,f=0.2,h=0.3,i=0.2,C1=0.0002,, C2=0.0005,C3=0.0003,A1=0.01,A2=0.02,andA3=0.03.Figure 17
Control profilesu1,u2,andu3, where a=0.1,b=0.1,c=0.1,e=0.3,f=0.2,h=0.3,i=0.2,C1=0.0002,C2=0.0005,C3=0.0003,A1=0.01, A2=0.02,andA3=0.03.
(a)(b)(c)


## 8. Conclusion

Our purpose in the present paper is to introduce insight into the Caputo fractional-order extension of the Lotka-Volterra model with an emphasis on immigration effects. Our main major objectives had been established which are as follows: The impact of adding small immigration terms to the prey or predators’ populations was deduced. The Grünwald–Letnikov method was applied to find numerical simulations for the formulated model. The effect of fractional-order derivative on the proposed model was studied. The stability problem was solved by adding the small immigration terms. Finally, the proposed small immigration terms can stabilize the natural populations of the one prey-two predators’ model. The fact that the findings forα=1 are in good agreement with the values chosen for α = 0.9, 0.85, and 0.75, which can be seen in all of the figures, demonstrating the effectiveness of the proposed technique in producing accurate results. All states have decreasing behavior as time marches towards a stable state. We conclude from these figures that the system is sensitive to the change in the fractional-order α and show how dependent it is on the history of the system. In the future, we will try to propose other new techniques for solving the stability problem concerned with the fractional-order one prey-two predators’ model. Also, we will consider cooperative phenomena; hence, a triple correlation in the equations will be taken into account. The equations would contain terms that depend simultaneously on all three variables.



---

*Source: 8795660-2023-02-22.xml*