Dynamics of a Predator–Prey Model with the Effect of Oscillation of Immigration of the Prey

In this article, the use of predator-dependent functional and numerical responses is proposed to form an autonomous predator–prey model. The dynamic behaviors of this model were analytically studied. The boundedness of the proposed model was proven; then, the Kolmogorov analysis was used for validating and identifying the coexistence and extinction conditions of the model. In addition, the local and global stability conditions of the model were determined. Moreover, a novel idea was introduced by adding the oscillation of the immigration of the prey into the model which forms a non-autonomous model. The numerically obtained results display that the dynamic behaviors of the model exhibit increasingly stable fluctuations and an increased likelihood of coexistence compared to the autonomous model.


Introduction
Organism population dynamics evolve through the changing sizes of populations. These dynamics take many different forms. Predator-prey interactions are primary forms that elucidate the population dynamics of many species. These interactions have received significant attention from many researchers [1][2][3][4][5][6]. Lotka [7] and Volterra [8] proposed the first model in this direction. Their model is considered to be the foundation stone for studying the dynamics of predator-prey interactions. Moreover, nonlinear differential equations form the mathematical basis for describing predator-prey interactions. The Lotka-Volterra model has been enhanced by many researchers to make it more realistic, and considering various ecological processes in doing so, see for instance [9,10].
Functional and numerical responses are considered as the central components for building predator-prey models; they also play a crucial role for describing the dynamic behaviors of these models [9]. The functional response is defined as the consumption rate of the prey by a predator, while the numerical response is the result of prey consumption on the predator density. In this context, the functional response was first suggested by Solomon [11] in 1949. A decade later, it was studied and classified into three types by the Canadian ecologist Holling, namely Holling types I, II and III [12]. The Holling type II was developed by Crowley and Martin [13], by adding terms that make the functional response depend on the prey and predator abundance, which takes into account the predator interference, so that makes it more realistic.
In the literature, few studies have investigated the Crowley-Martin type of responses in predator-prey models using different concepts to study the dynamics of these models. Upadhyay and Naji [14] used Holling type II and the Crowley-Martin type with a threespecies food chain model to discuss the dynamics of the model. Ali and Jazar [15] studied the global dynamics of a modified Leslie-Gower predator-prey model using the Crowley-Martin functional response and Holling type II as the numerical response. Stage-structured predator-prey models with the Crowley-Martin type have been widely considered [16][17][18] to discuss the dynamics of these models. Alebraheem and Abu-Hassan [19] investigated the seasonality in the predator-prey model with the Crowley-Martin type to identify complex dynamic behaviors. Mortoja et al. [20] used the Crowley-Martin type in a delayed predatorprey model involving disease in the prey population. Moreover, the Crowley-Martin type was used with a discrete predator-prey system to study the complex dynamics that lead to chaos [21,22]. In addition, the Crowley-Martin type has been investigated in stochastic predator-prey models [23,24] for studying the asymptotic properties of these models.
Different factors play an important role affecting the dynamic behaviors of predatorprey interactions, which may come from environmental or human factors such as refuge, migration, immigration, and oscillations. Immigration has effects on the stability of predator-prey interactions as shown in [25,26]. Oscillations that come from seasonal impacts play an important role affecting the dynamic behaviors of predator-prey interactions as shown in previous studies [27][28][29]. Despite the aforementioned efforts, no investigation has been conducted on oscillations in immigration that may affect the dynamics of predator-prey interactions.
In this article, Crowley-Martin-type functional and numerical responses are utilized to form an autonomous predator-prey model. The boundedness and validation of the model are presented. The dynamic behaviors of the model are studied by using some qualitative analysis, which involves stability analysis, and obtaining the coexistence and extinction conditions. In addition, a novel idea is introduced through adding the oscillation of the immigration of the prey into a predator-prey model, which forms a non-autonomous model. The comparison by using the numerical simulations between the two models are shown to present the change in dynamic behaviors and explaining the effects that come from investigating the oscillation of immigration of the prey.
The remainder of the article is set out as follows. Section 2 introduces the predatorprey model without the oscillation of immigration and proves the boundedness of the model. Section 3 presents the theoretical analysis of the model, which is divided to three subsections. Section 3.1 validates the model using Kolmogorov conditions. Section 3.2 introduces the existence of equilibrium points. Section 3.3 introduces the stability analysis of the model. Section 4 presents the numerical simulations that show the changing dynamic behaviors of this model with the oscillation of immigration of the prey. Section 5 presents a final discussion and conclusion.

Model and Boundedness without Oscillation of the Immigration
The autonomous predator-prey model with the Crowley-Martin functional and numerical responses is described as follows: which is subjected to the following initial conditions: Here, ζ represents the prey growth rate, k is the carrying capacity in the absence of predation, and µ is the predator death rate. The parameters of the functional and numerical responses are as follows: c denotes the efficiency of searching and capturing the predator, and γ indicates the efficiency of converting the consumed prey into predator births. Moreover, h represents the handling time, and σ is the magnitude of interference among predators. Theorem 1. For model (1) and any given positive initial value, the solution of t ≥ 0 will be in R 2 + and is ultimately bounded.
Proof. Suppose (x(t), y(t)) to be any solution of model (1) with initial conditions (2). Let The time derivative of function E is: Equation (3) becomes: Because all the parameters are positive, and the initiating solutions remain in R 2 + , the following could be assumed: Hence: Substituting (5) in (4), the following is obtained: which implies: Consequently, x(t) and y(t) are ultimately bounded.

Kolmogorov Analysis
Kolmogorov analysis was used to validate the predator-prey system and obtain the coexistence and extinction conditions. The Kolmogorov analysis is a group of conditions that are summarized and explained in Freedman [30] and Sigmund [31]. Therefore, the following conditions are used: (1) If the number of predators increases and the number of prey is fixed, then the prey and predator growth rate decrease. The factors that express these conditions are represented as follows: dN dy < 0.
Equation (8) can be interpreted as the predators competing for the same source. Therefore, these conditions are applied in model (1) as follows: The conditions (9) and (10) of model (1) are always negative because the values of the variables and parameters are positive.
(2) The model includes the environment carrying capacity. The following condition of model (1) is attained: (3) The model has a minimum prey value, even in the case of a small predator population.
The following condition of model (1) is attained: (4) The predators coexist with the prey if the following condition is satisfied: Thus, the following condition of model (1) is obtained: Condition (14) represents the coexistence condition. However, if this condition is not satisfied, then: kγc − kµh ≤ µ (15) and the predators become extinct. Many researchers have studied coexistence and extinction dynamics [32][33][34][35][36][37][38]. On account of their importance, the following definitions are presented: Definition 1. If x(0) > 0 and lim t→∞ x(t) > 0, then x(t) coexists. However, the geometrical definition is that each path of the differential equations is ultimately bounded afar from the coordinate planes.

Definition 2.
Extinction is analytically defined such that, if x(0) > 0 and lim t→∞ x(t) = 0, then x(t) become extinct. Meanwhile, the geometrical definition is that the path of the differential equations touches the coordinate planes.
From the Kolmogorov analysis, the following propositions are concluded: (14) is satisfied, then the predator coexists with the prey.

Existence of Equilibrium Points
Model (1) has three non-negative equilibrium points, which are biologically feasible because E 0 = (0, 0), E 1 = (k, 0), and E 2 = (x, y). Equilibrium points E 0 and E 1 exist without conditions. Equilibrium point E 2 exists as explained in Proposition 1. E 2 is obtained by finding the positive real root of the following algebraic equations: The following analysis is introduced to show the existence of x and y. From Equation (16), consider: The following is noticed: When x = 0, then y = ζ c−ζσ = y i . It is clear that y i > 0 under the following condition: (i) When y = 0, then x = k = x i . Then, it has x i > 0 always without any condition.
Then dy dx < 0, if it has the following condition: Using the above analysis, it was shown that isocline (18) passed through the points (x i , 0) and (0, y i ); and in Equation (18), y is a decreasing function of x under conditions (19) and (20).
From Equation (17), let: The following is noticed: (i) When y = 0, then x = µ γc−hµ = x ii . It is clear that x ii > 0 under the following condition: (21) is passing through the point (x ii , 0) under the condition (22), and it always has a positive slope, so in Equation (21), y increases as x increases.

Local and Global Stability Analysis
The qualitative analysis of differential equations is essential in population dynamics [39] because they follow the dynamic behaviors of populations and an exact solution is quite difficult [40]. The local and global stability of equilibrium points E 0 , E 1 , and E 2 are studied using the following theorems.
Proof. The Jacobian matrix of E 0 is yielded by Using the corresponding Jacobian matrix, it was obtained that, with eigenvalues λ 1 = ζ > 0, E 0 is always positive, and with λ 2 = −µ < 0, it is always negative. Therefore, E 0 is a saddle point.

Theorem 3.
Predator-free equilibrium point E 1 = (k, 0) is locally asymptotically stable and subjected to the following condition: Proof. The Jacobian matrix of E 1 is: Through the corresponding Jacobian matrix, it is attained that the eigenvalue λ 1 = −ζ and is always negative. Meanwhile, λ 2 = −µ + γck 1+hk will be negative if the condition (24) is satisfied.
Therefore, E 1 is locally asymptotically stable under condition (24). However, if the condition (24) is not satisfied, then E 1 is a saddle point.

Theorem 4.
Coexistence equilibrium point E 2 = (x, y) is locally asymptotically stable subjected to the following condition: Proof. The Jacobian matrix of E 2 is: Using the corresponding Jacobian matrix, E 2 is locally asymptotically stable under condition (25). (25) is violated, then the coexistence equilibrium point E 2 = (x, y) presents a stable fluctuated dynamic behavior.

Corollary 2. If condition
Theorem 5. Predator-free equilibrium point E 1 = (k, 0) is globally asymptotically stable in the interior of R 2 + under the condition (24).
Proof. Define the Lyapunov function of E 1 : where L is a continuously differentiable real valued function defined in R 2 + . Therefore, we have: If condition (24) is satisfied, then dL dt < 0 for any point in R 2 + .
Theorem 6. Coexistence equilibrium point E 2 is globally asymptotically stable in the interior of R 2 + under condition (25).
Proof. Consider that G(x, y) = 1 xy , where G represents a Dulac function. It is the C 1 function in the interior of R 2 and: T 2 (x, y) = −µy + γcxy 1 + hx + σy + hσxy . Thus: As x > 0, y > 0, and all the parameters are positive, then ∇ < 0 under condition (25). By the Bendixson-Dulac criterion, E 2 is globally asymptotically stable in the interior of R 2 + under the condition (25).

Model with Oscillation of Immigration
The non-autonomous predator-prey model with Crowley-Martin functional and numerical responses with oscillation of immigration of the prey is described as follows: which is subjected to the following initial conditions: where i the number of prey immigrants [25,26], represents the degree of fluctuation. The parameter w is the angular frequency of the fluctuations.

Numerical Simulations
The numerical simulations of models (1) and (32) were performed to show the change in dynamic behaviors and explain the effects that come from investigating the oscillation of the immigration of the prey. "NDsolve" command in the MATHEMATICA 11.3 software package was used to solve the models numerically for different sets of the hypothetical values of the parameters. They were assumed to satisfy the theoretical side of each case, but the values of the initial conditions were fixed for all cases. In Section 3, the theoretical analysis presents two different dynamics of the model (1), so two different sets of the hypothetical values of the parameters were selected for representing two different dynamic behaviors that were steady state and fluctuated, respectively. There is an important question that can be asked in this context: how does adding the oscillation of immigration of the prey affect the dynamic behaviors? Figures 1-10 illustrate the different kinds of graphs that are used to explain these cases; time series figures and zero-growth isoclines with phase plane trajectory figures are plotted to show the dynamic behaviors and trajectories comprehensively, as well as the cross-sectional picture of zero-growth isoclines with phase plane trajectories for each figure is introduced to give clear picture of the dynamic behavior.              The second set of parameter values is considered as follows:              The initial conditions for all cases are supposed as follows: x(0) = 3, y(0) = 2.
The first set of parameter values is considered as follows: By using the first set of parameter values without considering the oscillation of the immigration of the prey, all the Kolmogorov conditions are satisfied. Thus, according to Proposition 2, the unique positive coexistence equilibrium point is (3.39753, 3.07037). According to Theorem 4 and Theorem 6, it is locally and globally asymptotically stable since it satisfies condition (25), see Figure 1. Therefore, the dynamic behavior of the model (1) in the first set is steady state coexistence according to Proposition 2 and Theorem 6.
However, when investigating the oscillation of immigration of the prey, the dynamic behaviors of the model (32) tend to exhibit stable fluctuated and the fluctuations increase because of the increasing immigration parameter and the oscillation parameter that exists within the immigration, which is displayed through Figures 2-4. On the other side, when the value of the immigration parameter (i) increases, the joint equilibrium point of predator and prey increases and moves away from the axes where the predator and prey isoclines cross as shown through Figures 2b, 3b, 4b and 5b. It was interpreted that the likelihood of the coexistence of the predator-prey system increases as the value of i increases. Moreover, in the last case (i.e., i = 3.0), when neglecting the oscillation parameter i.e., = 0, the likelihood of the coexistence of the predator-prey system increases more because the joint equilibrium point of the predator and prey increases and moves away from the axes where the predator and prey isoclines cross as displayed through Figure 5b, in addition, the densities of the prey and predator increase as shown in Figure 5a.
Using the second set of parameter values without considering the oscillation in immigration, the dynamic behavior of model (1) is fluctuated coexistence, but close to axes, which increases the likelihood of extinction as depicted in Figure 6. This case corresponds with Proposition 1 and Corollary 2. In addition, all the Kolmogorov conditions are satisfied and the dynamic behavior stably fluctuates.
However, when investigating the oscillation of immigration, the dynamic behavior of the model (32) still fluctuated with the increasing immigration parameter and the oscillation parameter that exists within the immigration. Although the dynamic behavior is fluctuated, it moves away from the axes as shown through Figures 7b, 8b, 9b and 10b, and it is interpreted that the likelihood of coexistence of the predator-prey system increases as the value of i increases. Moreover, in the last case (i.e., i = 3.0), when neglecting the oscillation parameter, i.e., = 0, the likelihood of the coexistence of the system increases because the joint equilibrium point of the predator and prey increases and moves away from the axes where the predator and prey isoclines cross as displayed through Figure 10b, in addition, the densities of the prey and predator increase as shown in Figure 10a.
In Figures 9 and 10, it was noticed that the dynamic of the system tends to exhibit more stably and coexist as shown through the Figure 9b or Figure 10b that present the changing shape of the predator isocline, which can be interpreted that the dynamic of the system is in state of disequilibrium when using the second set of parameter values, but when adding the oscillation of immigration of the prey into the system, makes the system tend to exhibit a more stable coexistence.

Discussion and Conclusions
In this study, one of the most preferable and complicated predator-dependent functional and numerical responses is investigated to form an autonomous predator-prey model. The functional and numerical responses are namely Crowley-Martin type. The Crowley-Martin type is preferable because it was developed based on the Holling type II, which is well documented in the literature. The main dynamic behaviors of model (1) are studied. The boundedness of the model (1) is proved in Theorem 1.
The Kolmogorov analysis is employed to endorse the predator-prey model (1) and to identify the coexistence and extinction conditions. Furthermore, the local stability was analyzed by using the Jacobian matrix. The eigenvalues were determined to obtain the local stability conditions. Furthermore, the global stability is proved by using a Lyapunov function and the Bendixson-Dulac criterion under the same local stability conditions.
The theoretical results of this study are explained as follows. Our model has three non-negative equilibrium points, which is biologically feasible. However, the coexistence equilibrium point is the most important point because it represents the predator-prey interaction. This exists under certain conditions, as explained through Proposition 1. The stability conditions are obtained for each equilibrium point, as explained in Theorems 2-6. Corollary 2 explains that the coexistence equilibrium point has periodic dynamics. The coexistence and extinction conditions are identified, as shown in Proposition 2 and Corollary 1.
The oscillation of immigration of the prey is investigated to form a non-autonomous model. To the best of our knowledge, this is the first study to investigate the oscillation of immigration in predator-prey models; this makes the model more effective and realistic due to the existence of real factors such as environmental and human factors that may cause a fluctuation of immigration. The numerical simulations were performed by taking two different sets of the hypothetical values of the parameters to represent two different dynamics that are steady state and fluctuated, respectively. The numerical results showed that there are new remarkable outcomes that can be summarized as follows: When using the first set of parameter values, the following can be concluded: • Without consideration the oscillation of the immigration of the prey, the dynamic behavior of the system is in steady state coexistence. This coincides with the theoretical analysis of Proposition 2 and Theorem 6. • However, when investigating the oscillation of immigration of the prey, the dynamic behavior of the system tends to exhibit stable fluctuations which increase because of the increase in the immigration parameter and the oscillation parameter that exists within the immigration.

•
The likelihood of coexistence of the system increases as the value of immigration parameter increases. In addition to, when neglecting the oscillation parameter i.e., = 0, the likelihood of coexistence of the system increases more.
When using the second set of parameter values, the following is concluded: • Without consideration the oscillation of immigration of the prey, the dynamic behavior of the system is in fluctuated coexistence. This coincides with the theoretical analysis in Proposition 1 and the Corollary. • However, when investigating the oscillation of the immigration of the prey, the dynamic behavior of the system tends to exhibit stable fluctuations which increase because of the increase in the immigration parameter and the oscillation parameter that exists within the immigration.

•
The likelihood of the coexistence of the system increases as the value of the immigration parameter increases. In addition, when neglecting the oscillation parameter, i.e., = 0, the likelihood of coexistence of the system increases more.

•
The dynamic of the system tends to exhibit a more stable coexistence as the immigration parameter increases.
Some of the results of this study are consistent with previous findings on the effect of small immigration in a predator-prey system, where it demonstrates that the system stabilizes with adding positive immigration [26], as in the latter cases of the first and second sets in this study. However, the dynamic behaviors exhibit increasing fluctuations, but they display stable fluctuations and the coexistence of the system was paid more attention in this study, where the likelihood of the coexistence of the system increases as the value of the immigration parameter increases. In nature, there are some ecological interpretations may support the obtained outcomes in this study such as the rescue effect [41], which is a phenomenon explaining the increase in the likelihood of coexistence through immigration [42,43].