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

,


Introduction
Ecology is the scientifc study that relies on the distribution and abundance of living things and how these factors are afected 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 diferent felds 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 fnding the analytical and numerical solutions of such biomathematical models, we can get a clear understanding of the underlying phenomena.Ten, the control strategies can be designed.
Te interaction between the prey and predators is the main feature that reveals the behavior of the prey and predators.Tere 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].
Te dynamical behavior of an exploited system consisting of two preys and a predator, that is, being harvested, is investigated in [11].Te 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].Te authors in [13] examined the predator-preypathogen relationship, in which the predator afects the infection's transmission rate in its prey.A predator-preypathogen model in which the predator is a natural specialist and infected prey can enter refugia of a fxed size to avoid predator attack was proposed in [14].Te bifurcation analysis of a predator-prey model with prey refuge is discussed in [15].Recently, the classical Lotka-Volterra (LV) model has been modifed by several researchers by studying several factors such as disease, Allee efects, and fear mechanisms.For the classical LV model, the mathematical properties indicate either divergent extinction or cyclic oscillation of one species, see [16].Te 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.Te 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 efect.Several felds, such as biology, fnance, economics, engineering, and medicine, have used fractional calculus to solve their problems; see for example [19] and the references therein.Tere exist several defnitions for fractional-order derivatives.Te most popular defnition of fractional-order derivative is the Caputotype 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 diferential equations with exogenous input and proposed the generalized Mittag-Lefer input stability conditions.Rashid et al. [26] proved a multiple fxed-point result for partially ordered s-distance spaces under (θ, ϕ, ψ)-type weak contractive condition.In [27], a fractionalorder 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.Te 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-Lefer 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][33][34][35].Tahara et al. [36] studied the efect of immigration on a Lotka-Volterra model.However, no work has been performed on the Caputo fractional-order Lotka-Volterra model with the immigration efect.Our proposed model is an extension to the model in [36], where the proposed model is a one prey-two predators' model.
Te 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 efect in this paper.Ten, we study the efect of adding a small immigration term to the prey or predator populations.Te extension is both the inclusion of a second predator and the use of Caputo derivatives.Te reason for selecting this particular Caputo defnition is that, unlike the Riemann-Liouville fractional derivative, the Caputo fractional derivative allows for the inclusion of initial conditions in the model.Te Caputo fractional derivative of a constant is also zero, in contrast to the Riemann-Liouville fractional derivative.Te 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 modifcations.Tis study comes to the conclusion that very little immigration into either the prey population or the predators' population stabilizes the Lotka-Volterra system.All fuctuations 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.Te following few points highlight the proposed research's novelty: (1) A newly developed fractional model, defned in the Caputo sense, for examining the new dynamics of the one prey-two predators' model with the efect of immigration.
(3) To validate the theoretical facets and results of the suggested model, an efective numerical technique is adopted.(4) Te obtained results demonstrate the efectiveness of the suggested model in providing some new insights into the dynamical behavior of the proposed model.(5) An optimal control strategy is discussed.Tis paper is divided as follows: Section 2 presents the fundamentals and basic defnitions that will be utilized later in other sections.Section 3 displays the fractionalorder Lotka-Volterra model with immigration efect.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 2

Journal of Mathematics
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.

Fundamentals and Definitions
Some of the fundamental defnitions that will be utilized later in other sections are provided in this section, see [37,38].
Defnition 1. Te Riemann-Liouville fractional integral for a continuous function f: [0, +∞) ⟶ R is given by the following equation: ( Defnition 2. Te Riemann-Liouville fractional derivative for a continuous function f: [0, +∞) ⟶ R can be represented by the following equation: Defnition 3. Te Caputo fractional derivative for a continuous function f: [0, +∞) ⟶ R is defned by the following equation: Grünwald-Letnikov (GL) defnition of fractional derivatives is given by the following equation: where [.] refers to the integer part.By applying relation (5) formulated from relation (4), we can calculate the fractional-order derivatives numerically.Te numerical approximation of α-th derivative at the points kh, (k �1, 2, ...) can be written as follows [39]: where L m refers to the length of memory, t k � kh, h is step size, and c (α) i (i � 0, 1, . .., k) are binomial coefcients.Te following expression can be used to calculate c (α) i as follows: Te nonlinear fractional diferential equation in the form a D α t f(t) � v(x(t), t) has a general numerical solution that can be expressed by the relations (5) and (6) as follows: A short memory principle can be used.Ten, the lower index of the sums in relations (7) will be p � 1 for k < (L m /h) and p � k − (L m /h) for k > (L m /h).Scherer et al. [40] studied the existence of bounded unique solutions for the Grünwald-Letnikov defnition 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.

Proposed Fractional Model
Here, we propose the Caputo fractional Lotka-Volterra model with few prey and predators immigrants to refect the efect of the memory.Also, we study the efect of adding an immigration term D(x) to the prey population or adding immigration terms G(y) and J(z) to the predators populations.
Te Caputo fractional one prey-two predators model with few immigrants can be represented as follows: where a, b, c, e, f, h, and i are positive real constants.Te parameter a refers to the reproduction rate of the prey, b refers to the rate at which the frst predator consumes the prey, and c refers to the rate at which the second predator consumes the prey.Te parameter e refers to the birth rate of the frst predator for each prey captured, and f refers to the mortality rate of the frst predator.Te 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.Te prey population is represented by x, the frst predator population is represented by y, and the second predator population is represented by z.C 0 D α t refers to the Caputo fractional derivative of order α.Te proposed model is a generalization to the model introduced in [36].Te immigration function can be represented in three ways as follows: where d represents the number of prey immigrants, g represents the number of frst 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 frst predator immigrants, and j/z represents the proportion of second predator immigrants.We proposed the following assumptions for equations ( 9)-( 11): x ≠ 0 when D(x) � d/x, y ≠ 0 when G(y) � g/y, and z ≠ 0 when J(z) � j/z.In this paper, to study the efect of the immigration factor on the Caputo fractional one prey-two predators model (system (8)), we consider the following six cases:

Existence, Uniqueness, and Non-Negativity
We established several fndings regarding the existence and uniqueness of solutions for system (8) using fxed point theory fndings, such as the Schauder and Banach fxed point theorems.We follow [43] in this section.Te following compensation is used for the right-hand sides of system (8): Tis yields the following system: t on both sides of equation ( 13) and using the initial conditions, one has the following equation: 4 Journal of Mathematics We consider the following assumptions to establish the results:

Tere are real positive values
(ii) Tere are constants E ψ 1 , E ψ 2 , E ψ 3 > 0, such that for (x, y, z), (x, y, z) ∈ B, we have the following equation: Theorem 1. Te considered system (8) has at least a solution under assumption (i).
we have the following equation: Journal of Mathematics and similarily Tis implies that Ξ(V) ⊂ V. Also ψ 1 , ψ 2 , ψ 3 are continuous so is Ξ.For t 1 < t 2 ∈ [0, τ], we have the following equation: It is clear that the right side of equation ( 20) goes to zero as t 2 ⟶ t 1 and similarly we obtain the following equation: Consequently, Terefore, 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 fxed point, the operator Ξ has at least a fxed point.Hence, system (8) has at least a solution.
and the assumption (ii) is verifed, then system (8) has a unique solution.
Proof.Suppose (x, y, z), (x, y, z) ∈ B, then consider 6 Journal of Mathematics Ten, from equation ( 22), we obtain Hence Ξ is a contraction operator and the system under study (8) has a unique solution.
Only non-negative solutions are of importance to us because of the biological signifcance.Te following theorem proves that the solutions of model ( 8) are non-negative.Suppose R + is the set of all non-negative real numbers and Theorem .For the model ( 8), all the solutions that start in Λ + are non-negative.
Proof.Firstly, we prove that the solutions x(t) which start in Λ + are non-negative, that is, x(t) ≥ 0 for all t ≥ t 0 .Suppose that is not true, then there exists a constant t 1 ≥ t 0 such that Based on equation ( 23) and the frst equation of system (8), we have the following equation: According to Lemma 1, we have x(t + 1 ) ≥ 0, which contradicts the fact x(t + 1 ) < 0. Terefore, we have x(t) ≥ 0 for all t ≥ t 0 .Using the same manner as above, we have y(t) ≥ 0, z(t) ≥ 0 for all t ≥ t 0 .

Stability Analysis
We can fnd the equilibrium points for the proposed model ( 8) by solving the following system: Te formulation of the previous system ( 25) is based on the values of the functions D(x) and G(y), and Tables 1-4 show the equilibrium points and the corresponding eigen values for all the proposed cases of D(x), G(y), and J(z).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 Tables 1-4, Eig(J) 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 Eig(J) 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: As stated above, the equilibrium points of system (8) can be found by solving the following equation: f(x) � 0. Tese points are locally asymptotically stable if all eigenvalues λ i of the Jacobian matrix J � zf/zx evaluated at the equilibrium points satisfy the Matignon conditions: Te stable and unstable regions for 0 < α < 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.

Analysis of Results.
Trough this paper, we substitute x(0) � 5, y(0) � 5, and z(0) � 10.A time simulation of the proposed model is presented by applying the Grünwald-Letnikov method.First, we show the efect of changing the fractional order on the one prey-two predator model without immigration (original case) (i.e., D(x) � 0, G(y) � 0, J(z) � 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 afected by changing the value of the fractional-order α.Te phase space of the one prey-two predators system is shown in Figure 3.
For Case A1, we have D(x) � d, G(y) � 0, and 4 depicts how the population dynamics are afected by changing the value of the fractional-order α.Te phase space of the one prey-two predators is shown in Figure 5.
For Case C1, we have D(x) � d, G(y) � 0, and J(z) � 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 afected by changing the value of the fractional-order α.Te phase space of the one preytwo predators system is shown in Figure 7.
For the Case E1, we have D(x) � 0, G(y) � g/ y, and J(z) � 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 afected by changing the value of the fractional-order α.Te phase space of the prey-two predators system is shown in Figure 9.
For the Case B2, we have D(x) � 0, G(y) � − g, and J(z) � 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 afected by changing the value of the fractional-order α.Te phase space of the one prey-two predators system is shown in Figure 11.
For the Case F2, we have D(x) � 0, G(y) � 0, and J(z) � − 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 afected by changing the value of the fractional-order α.Te phase space of the one prey-two predators system is shown in Figure 13.8 Journal of Mathematics Journal of Mathematics Journal of Mathematics and α � 0.75, 0.85, 0.9, 1; and 3D-plot of x(t), y(t), and z(t) for α � 1.    16 Journal of Mathematics systems.Keep in mind that if the prey goes extinct, the predator population will also go extinct.Tis does not, however, imply that if the predator population goes extinct, the prey will follow suit.the population of prey will continue to exist over time and even increase.Tere 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.Te current research demonstrates that extremely few migrants, whether they are predators or prey, produce stability in the LV systems.Te biological meaning of immigration changes when d (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 fuctuations in the population of predators and prey will be averaged out by the addition of positive immigration.Additionally, the fndings show that a positive immigration component is sufcient 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.Tese few migrant species are enough to keep the prey-predator models stable.Tis 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 α.Te observed behavior from the numerical simulations shows that the states have a decreasing behavior as time is increased while the fractional-order α is decreased.Te efect of the fractional derivative on the behavior of the dynamics of the proposed model is noticed to show its stabilization efect.

Optimal Control Strategy
Here, we proposed three time-dependent control techniques in the same time.We calculated the optimal values for u 1 , u 2 , and u 3 that would reduce the oscillations in (x, y, z).Additionally, the costs are minimized.For more illustrations for diferent optimal control strategies, see [47][48][49] and the references therein.Te proposed system for the optimal control is as For U � (u 1 (t), u 2 (t), u 3 (t)) ∈ U, where the objective functional is defned by the following equation: where Here, A 1 represents the weight for the number of the frst predator, A 2 represents the weight for the number of the second predator, A 3 represents the weight for the number of the prey, and C i (i � 1, 2, 3) stand for weight coefcients and express the unit costs of the controls at diferent levels.We will determine the required conditions for the optimal control U * based on Pontryagin's maximum principle [50].Theorem 4. Let (x, y, z) be the state solutions of problem ( 27) and ( 28) and (u * 1 , u * 2 , u * 3 ) be the optimal controls, then there are adjoint functions satisfying the following equation: Journal of Mathematics  with the transversality conditions λ i (T) � 0 and i � 1, 2, 3. Furthermore, the optimal controls (u * 1 , u * 2 , u * 3 ) are as follows: Proof.Using Pontryagin's maximum principle, if U * ∈ U is considered an optimal solution ( 27) and ( 28), then there exists a continuous mapping λ: [0, T] ⟶ R, where λ(t) � (λ 1 (t), λ 2 (t), λ 3 (t)) is the adjoint vector, which is obtained from the following equation:   Te Hamiltonian function ϕ is as follows: Te adjoint system can be obtained as a result.To obtain the characterization of the optimal control U * , we use the following equation: 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.Figure 14 shows the efect of applying the control signals u 1 , u 2 , and u 3 on the state x(t).Te solid lines show the state x(t) without applying the control signals u 1 , u 2 , and u 3 .Te dashed lines show the state x(t) after applying the control signals u 1 , u 2 , and u 3 .Figure 15 shows the efect of applying the control signals u 1 , u 2 , and u 3 on the state y(t).Te solid lines show the state y(t) without applying the control signals u 1 , u 2 , and u 3 .Te dashed lines show the state z(t) after applying the control signals u 1 , u 2 , and u 3 .Figure 16 shows the efect of applying the control signals u 1 , u 2 , and u 3 on the state z(t).Te solid lines show the state z(t) without applying the control signals u 1 , u 2 , and u 3 .Te dashed lines show the state z(t) after applying the control signals u 1 , u 2 , and u 3 .It is clear that applying the control signals reduces the oscillations.Figure 17 shows the control profle u 1 , u 2 , and u 3 at diferent values of α � 1, 0.9, and 0.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 efects.Our main major objectives had been established which are as follows: Te impact of adding small immigration terms to the prey or predators' populations was deduced.Te Grünwald-Letnikov method was applied to fnd numerical simulations for the formulated model.Te efect of fractional-order derivative on the proposed model was studied.Te 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.Te fact that the fndings 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 fgures, demonstrating the efectiveness of the proposed technique in producing accurate results.All states have decreasing behavior as time marches towards a stable state.We conclude from these fgures 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.Te equations would contain terms that depend simultaneously on all three variables.

Table 1 :
Equilibrium points and eigen values for the original case, case A1, case B1, and case C1.

Table 2 :
Equilibrium points and eigen values for case D1, case E1, and case F1.

Table 3 :
Equilibrium points and eigen values for case A2, case B2, and case C2.

Table 4 :
Equilibrium points and eigen values for case D2, case E2, and case F2.