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.

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 [3235]. 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: Section 2 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 function is given by the following equation:

Definition 2. The Riemann–Liouville fractional derivative for a continuous function can be represented by the following equation:

Definition 3. The Caputo fractional derivative for a continuous function is defined by the following equation:Grünwald–Letnikov (GL) definition 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. The numerical approximation of -th derivative at the points ( = 1, 2, ...) can be written as follows [39]:where refers to the length of memory, is step size, and (i = 0, 1, …, are binomial coefficients. The following expression can be used to calculate as follows:The nonlinear fractional differential equation in the form has a general numerical solution that can be expressed by the relations (5) and (6) as follows:A short memory principle can be used. Then, the lower index of the sums in relations (7) will be p = 1 for and for . 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 term to the prey population or adding immigration terms to the predators populations.

The Caputo fractional one prey-two predators model with few immigrants can be represented as follows:where and are positive real constants. The parameter refers to the reproduction rate of the prey, refers to the rate at which the first predator consumes the prey, and refers to the rate at which the second predator consumes the prey. The parameter refers to the birth rate of the first predator for each prey captured, and refers to the mortality rate of the first predator. The parameter refers to the birth rate of the second predator for each prey captured, and refers to the mortality rate of the second predator. The prey population is represented by the first predator population is represented by and the second predator population is represented by 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:where represents the number of prey immigrants, represents the number of first predator immigrants, and represents the number of second predator immigrants. Also, represents the proportion of prey immigrants, represents the proportion of first predator immigrants, and represents the proportion of second predator immigrants. We proposed the following assumptions for equations (9)–(11): when , when , and when .

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 ()Case B1: immigrants of the first predator ()Case C1: immigrants of the second predator ()Case D1: few immigrants of the prey ()Case E1: few immigrants of the first predator ()Case F1: few immigrants of the second predator ()

We also handle the following six cases of the immigration factor:Case A2: migrants of the prey (i.e., )Case B2: migrants of the first predator (i.e., )Case C2: migrants of the second predator (i.e., )Case D2: few migrants of the prey (i.e., )Case E2: few migrants of the first predator (i.e., )Case F2: few migrants of the second predator (i.e., )

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):

This yields the following system:

Clearly, upon using integral operator on both sides of equation (13) and using the initial conditions, one has the following equation:

Let be Banach space equipped with the norm . We defined an operator by the following equation:

We consider the following assumptions to establish the results:There are real positive values , such that(ii)There are constants , such that for , we have the following equation:

Theorem 1. The considered system (8) has at least a solution under assumption (i).

Proof. Suppose is a convex and closed set, then to prove that is continuous, suppose , we have the following equation:and similarilyThis implies that . Also are continuous so is For , we have the following equation:It is clear that the right side of equation (20) goes to zero as and similarly we obtain the following equation:Consequently,
as . Hence, we have as .
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 and the assumption (ii) is verified, then system (8) has a unique solution.

Proof. Suppose , then considerThen, from equation (22), we obtain Hence is a contraction operator and the system under study (8) has a unique solution.

Lemma 1 (See [44]). Let Suppose that and . If , then is a nondecreasing function for each . If , then is a nonincreasing function for each .

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 .

Theorem 3. For the model (8), all the solutions that start in are non-negative.

Proof. Firstly, we prove that the solutions which start in are non-negative, that is, for all . Suppose that is not true, then there exists a constant such thatBased on equation (23) and the first equation of system (8), we have the following equation:According to Lemma 1, we have , which contradicts the fact . Therefore, we have for all . Using the same manner as above, we have for all .

5. Stability Analysis

We can find the equilibrium points for the proposed model (8) by solving the following system:

The formulation of the previous system (25) is based on the values of the functions and , and Tables 14 show the equilibrium points and the corresponding eigen values for all the proposed cases of , , and . 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 14, refers to the eigen values of the Jacobian matrix that can easily be computed by solving for We have used the symbol 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: . These points are locally asymptotically stable if all eigenvalues of the Jacobian matrix evaluated at the equilibrium points satisfy the Matignon conditions:

The stable and unstable regions for are shown in Figure 1. It is clear that there are many eigen values in Tables 14. Hence, in what follows and based on Lemma 2, we will discuss the stability of some equilibrium points in Tables 14.(i)From Table 1, the corresponding eigen values for the trivial equilibrium point are . Hence, for all for all and for all Consequently, the trivial equilibrium point is unstable.(ii)From Table 1, the corresponding eigen values for the equilibrium point are . If then for all , , and Consequently, the equilibrium point is unstable.(iii)From Table 1, the corresponding eigen values for the equilibrium point are If , then for all , and for all . Consequently, the equilibrium point is asymptotically stable.(iv)From Table 4, the corresponding eigen values for the equilibrium point are . Hence, for all , and . Consequently, the equilibrium point 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.

6. Numerical Simulation

6.1. Analysis of Results

Through this paper, we substitute . 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., ) By taking and 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.

For Case A1, we have By taking and 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.

For Case C1, we have By taking and 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.

For the Case E1, we have By taking and 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.

For the Case B2, we have By taking and 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.

For the Case F2, we have By taking and 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.

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 when (or or ) and (or or ) 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 213. 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 213 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 for that would reduce the oscillations in . Additionally, the costs are minimized. For more illustrations for different optimal control strategies, see [4749] and the references therein. The proposed system for the optimal control is as

For wherethe objective functional is defined by the following equation:where

Here, represents the weight for the number of the first predator, represents the weight for the number of the second predator, represents the weight for the number of the prey, and 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 based on Pontryagin’s maximum principle [50].

Theorem 4. Let be the state solutions of problem (27) and (28) and be the optimal controls, then there are adjoint functions satisfying the following equation:with the transversality conditions Furthermore, the optimal controls are as follows:

Proof. Using Pontryagin’s maximum principle, if is considered an optimal solution (27) and (28), then there exists a continuous mapping where is the adjoint vector, which is obtained from the following equation:that is,The Hamiltonian function is as follows:The adjoint system can be obtained as a result. To obtain the characterization of the optimal control , 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 effect of applying the control signals on the state The solid lines show the state without applying the control signals The dashed lines show the state after applying the control signals Figure 15 shows the effect of applying the control signals on the state The solid lines show the state without applying the control signals The dashed lines show the state after applying the control signals Figure 16 shows the effect of applying the control signals on the state The solid lines show the state without applying the control signals The dashed lines show the state after applying the control signals It is clear that applying the control signals reduces the oscillations. Figure 17 shows the control profile at different values of and 0.8.

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 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.

Data Availability

The data used to support the findings of this study are available from the corresponding author on request.

Conflicts of Interest

The authors declare that they have no conflicts of interest.