Global dynamics of a predator-prey system with immigration in both species

Electronic Research Archive http://www.aimspress.com/journal/era ERA, 32(2): 762–778. DOI: 10.3934/era.2024036 Received: 13 November 2023 Revised: 26 December 2023 Accepted: 03 January 2024 Published: 10 January 2024 Research article Global dynamics of a predator-prey system with immigration in both species Érika Diz-Pita* Departamento de Estatística, Análise Matemática e Optimización, Universidade de Santiago de Compostela *Correspondence: Email: [email protected]. Abstract: In nature, the vast majority of species live in ecosystems that are not isolated, and the same is true for predator-prey ecological systems. With this work, we extend a predator-prey model by considering the inclusion of an immigration term in both species. From a biological point of view, that allows us to achieve a more realistic model. We consider a system with a Holling type I functional response and study its global dynamics, which allows to not only determine the behavior in a region of the plane R2, but also to control the orbits that either go or come to infinity. First, we study the local dynamics of the system, by analyzing the singular points and their stability, as well as the possible behavior of the limit cycles when they exist. By using the Poincaré compactification, we determine the global dynamics by studying the global phase portraits in the positive quadrant of the Poincaré disk, which is the region where the system is of interest from a biological point of view. Keywords: predator-prey; immigration; stability; global dynamics; phase portrait 1. Introduction and statement of the main result In the field of biosystems, dynamical systems and differential equations have been widely used for decades. Numerous works have attempted to solve various problems in ecology, biology, and medicine, from the classic Lotka-Volterra or SIR models to some very recent works such as [1], where the authors studied vegetation patterns with the aim of either preserving or restoring them using a vegetation model with non-local root interaction, or [2], where the authors relate human activities to the evolution of ecosystems, trying to provide solutions for biodiversity conservation. Since the beginnings of biomathematics, predator-prey systems have been a special area of great interest to researchers. Advances in these models are continuously being made; some recent studies are highlighted in [3–10]. Although there are many contributions in the literature on predator-prey models, due to the complexity of the real systems being represented, there are many features that still need to be studied in greater detail. For example, it is important to consider the effects of the presence 763 of some number of immigrants, as most systems in nature are not isolated; moreover, in all major branches of the animal kingdom, there are numerous migrating species, from crustaceans or insects to large mammals [11]. In models of very different scales, this migratory phenomenon must be taken into account , since it affects animals of various sizes, from small zooplankton species (1 mm in length) to large blue whales (up to 27 meters in length) [12,13]. In recent years, numerous studies have analyzed the effects of the presence of immigrants on different populations and species, under different hypotheses and with different mathematical tools and models.There are some works that used delay equations because delayed migration can occur when the individuals encounter some barriers, see for example [14–17]. Other recent works considered fractional order models, as in [18]. Among these works, some have been selected, analyzed and compared in [19]. We highlight the work of Sugie [20], in which the authors included the effect of a constant immigration rate that affected the prey species in the classical Rosenzweig-MacArthur model. This allowed them to state some ecological conclusions that highlighted the fact that the inclusion of immigration allows for the coexistence of species, and that considering immigration can be important and have considerable effects even in simple models. Another very recent work, is that of of Priyanka et al. [21], which considered immigration in both the prey and the predator. However, in that case, immigration is not constant but proportional to the population; additionally, immigration was combined with other characteristics such as harvesting, which was represented by a term with the same form to that of immigration. The work of Tahara et al. [22], analyzed different predator-prey models with immigration. In addition to different types of functional responses, they considered two representations of immigration: one as a constant and another as a function of the type c/x. Several numerical simulations were performed to show how the inclusion of immigration allowed for the stabilization of populations. Although there were some cases in which immigration in only one of the species was studied, no analytical study has studied the case with immigration in both species. Other works dealt with predator-prey models in higher dimensions. This is the case for the work of Mukherjee [23], in which they considered a predator prey model in dimension three thereby containing a predator species, a prey species, and a competitor of the prey. In this case, the considered immigration rate was also constant, but is only taken into account in the prey species. The author studied the existence and stability of the equilibria and showed that a Hopf bifurcation can occur under certain conditions. Additionally, we found some discrete-time systems that modeled this type of predator-prey ecosystem with immigration. As previously described, in these cases, they are limited to consider immigration in only one of the species. For example, in [24], the authors extended the study of the Holling type I discrete system considered in [25] by adding a constant prey immigration rate. They studied the equilibria and their topological classification, and showed the changes in the dynamics could be obtained with the inclusion of the immigration term. Moreover, bifurcations were also analyzed, and bifurcation diagrams together with phase portraits were included; however, they were obtained by numerical simulations and not directly from the analytical results. Motivated by these works in the literature, we sought to carry out a complete study of the global dynamics of a predator-prey system with immigration in both species. We will take the immigration rate as a constant, as it is done in many of the aforementioned works, which is reasonable from a biological Electronic Research Archive Volume 32, Issue 2, 762–778. 764 point of view. Then, we will be advance in the study of the dynamics in two ways: by adding a new characteristic (i.e., immigration) and by making a study of the dynamics through a compactification that is not limited to local behavior, thus providing us with the ability to understand the behavior at infinity. This has been done recently with other classical predator-prey models, as in [26]. We consider that this work can be very useful to study ecosystems in which both prey and predators present migratory behaviors, which is something that occurs in the nature. The mathematical study of the undermentioned system will allow one to predict and interpret the behavior of the species. On one hand, we would like to point out that we have not found any other work in which immigration is considered for both species, and this is one of the innovative aspects of the work presented. On the other hand, in general, the study of the global dynamics of the systems is not carried out, though this is something that allows us to have complete control of the dynamics, without limiting it to local behaviors; futhermore, this allows us to work in a compact region, and thus all the possible dynamic behaviors of a system can be classified. Therefore, we hope that this model, as well as the classification of its global dynamics, serves as a direct application to ecological problems, and as a starting point for future models that consider immigration in both the prey and predator populations In general, we can consider the following system: ˙x=rx −ax1+αy 1+hx1+α+c1, ˙y=bx1+αy 1+hx1+α−ny +c2, (1.1) where rrepresents the growth rate of prey, nis the death rate of predator, ais the rate of predation, b is the conversion rate of eaten prey into new predators, hand αare the functional response coefficients (which involve for example the handling time), and c1and c2are the immigration rates of prey and predator, respectively. All these real parameters are positive due to their biological meaning. To begin with, in the present work, we deal with the case where α=h=0, which corresponds with a Holling type I functional response. Our main result is as follows. Theorem 1.1. The global phase portrait of system (1.1) where α=h=0, in the positive quadrant of the Poincaré disc, is one of the following: Figure 1. Global phase portraits of system (1.1) with α=h=0 in the positive quadrant of the Poincaré disk. Electronic Research Archive Volume 32, Issue 2, 762–778. 765 This result is important as it allows us to determine that for any value of the parameters, the behavior will follow one of these two schemes. Furthermore, we emphasize that the given phase portraits are in the positive quadrant of the Poincaré disk (i.e., they are global phase portraits), which allows us to control the behavior near the infinity, and to determine how the orbits behave, either coming from or going to infinity. From an ecological point of view, in the first phase portrait, there is a singular point which represents the coexistence of both species, and we observe that regardless of the initial condition considered ( i.e., the initial number of prey and predators), the solution predicts an evolution to a singular point of coexistence. In the second phase portrait, there are two regions: the one inside the limit cycle and the one outside it. Within the limit cycle, the behavior is similar to that mentioned in the first phase portait, since given any initial condition in that region, the number of prey and predators will tend to the value at the equilibrium point. For the initial conditions outside the limit cycle, the solutions increasingly approach the limit cycle on the outside, which means that, in practice, there will be an oscillation in the number of prey and predators, with cyclic type of behavior. 2. Proof of the result 2.1. Finite singular points System (1.1), where h=α=0, corresponds to a Holling type I functional response: ˙x=rx −axy +c1, ˙y=bxy −ny +c2. (2.1) In order to make the reading more clear, we introduce the following notation: R=p4bc1nr +(bc1+ac2−nr)2. There are two general solutions for ˙x=˙y=0, with the following expressions −bc1−ac2+nr +R 2br ,bc1+ac2+nr +R 2an !and −bc1−ac2+nr −R 2br ,bc1+ac2+nr −R 2an !, which are always well defined since the parameters are positive and Ris a positive real number for all the values of the parameters. In any case, given the motivation behind the formulation of the system, we are only interested in the non-negative singular points; therefore, we analyze the sign and position of the points given by these expressions. The following results allow us to characterize the location of these singular points. Proposition 2.1. The region C =n(x,y)∈R2|x≥0,y≥0ois positively invariant. Proof. The fact that region Cis positively invariant can be deduced from the fact that ˙x|x=0=c1>0,and ˙y|y=0=c2>0,(2.2) Electronic Research Archive Volume 32, Issue 2, 762–778. 766 which means that the orbits of the system enter the region at any point of the boundary of C, when t moves on in the positive sense. □ Corollary 2.2. There are no singular points representing the survival of only one of the species. Proof. As the positive axes n(0,y)∈R2|y≥0oand n(x,0) ∈R2|x≥0oare not invariant lines, and the flow is transversal to them, the singular points of system (2.1) cannot be over the axes. □ Theorem 2.3. System (2.1) has exactly one positive singular point P= −bc1−ac2+nr +R 2br ,bc1+ac2+nr +R 2an !, for any values of the parameters. Proof. Let us prove that this point is always positive. For the second component, this is trivial, as all parameters are positive. For the first component, if we suppose that it is negative, then, −bc1−ac2+ nr +R<0. Then, R=p4bc1nr +(bc1+ac2−nr)2<bc1+ac2−nr; by squaring the expressions, 4bc1nr +(bc1+ac2−nr)2<(bc1+ac2−nr)2⇒4bc1nr <0, which contradicts the fact that all parameters are positive. An analogous reasoning allows us to prove that the other solution of ˙x=˙y=0 is never positive. If we suppose that it is positive, then −bc1−ac2+ nr −R>0; by squaring the expressions, one can obtain the following: (bc1+ac2−nr)2>4bc1nr +(bc1+ac2−nr)2⇒4bc1nr <0, which is a contradiction. Then, the system has always exactly one positive singular point. □ Proposition 2.4. The singular point P of system (2.1) is asymptotically stable and has the following phase portraits, depending on the parameters: 1) It is a stable focus if ((n+r)(bc1+ac2)+(n−r)(nr +R))2<16n2r2R. 2) It is a stable node if ((n+r)(bc1+ac2)+(n−r)(nr +R))2>16n2r2R. Proof. The Jacobian matrix of system (2.1) at the singular point Pis M= r−bc1+ac2+nr +R 2n a(bc1+ac2−nr −R) 2br b(bc1+ac2+nr +R) 2an −bc1+ac2+nr −R 2r  , and has the following eigenvalues: λ1=−1 4nr(ac2+bc1)(n+r)+(nr −R)(n−r) +q−16n2r2R+((n+r)(bc1+ac2)+(n−r)(nr −R))2!, Electronic Research Archive Volume 32, Issue 2, 762–778. 767 λ2=−1 4nr(ac2+bc1)(n+r)+(nr −R)(n−r) −q−16n2r2R+((n+r)(bc1+ac2)+(n−r)(nr −R))2. Then, the trace of the Jacobian matrix is as follows: tr(M)=r−bc1+ac2+nr +R 2n−bc1+ac2+nr −R 2r where the expression is always negative, and can be proved. First, r−bc1+ac2+nr +R 2n=−bc1+ac2−nr +R 2n is negative as bc1+ac2−nr +R=bc1+ac2−nr +p(bc1+ac2−nr)2+4bc1nr >0. Furthermore, taking into account that 4bc1nr +(bc1+ac2−nr)2=4bc1nr +(bc1−nr)2+a2c2 2+2(bc1−nr)ac2 =4bc1nr +b2c2 1+n2r2−2bc1nr +a2c2 2+2(bc1−nr)ac2 =b2c2 1+n2r2+2bc1nr +a2c2 2+2(bc1−nr)ac2 =(bc1+nr)2+a2c2 2+2(bc1+nr)ac2−4nrac2 =(bc1+ac2+nr)2−4nrac2, it holds that −bc1+ac2+nr −p(bc1+ac2−nr)2+4bc1nr 2r=−bc1+ac2+nr −p(bc1+ac2+nr)2−4bc1nr 2r which is negative as bc1+ac2+nr >p(bc1+ac2+nr)2−4nrac2.Then, we have that the trace is always negative for any values of the parameters. The determinant of the Jacobian matrix is as follows: λ1λ2=qb2c2 1+(ac2−nr)2+2bc1(ac2+nr), which is positive. If the radicand in the expressions of λ1, λ2is negative ( i.e., if ((n+r)(bc1+ac2)+(n−r)(nr +R))2<16n2r2R), then the two eigenvalues are complex, and the singular point Pis a stable focus as λ1+λ2=2Re(λ1)=2Re(λ2)<0. Then, Pis asymptotically stable. If the radicand is positive ( i.e., if ((n+r)(bc1+ac2)+(n−r)(nr +R))2>16n2r2R), then the eigenvalues are real; and since λ1λ2>0, the singular point is a node, and taking into account that λ1+λ2<0, then both eigenvalues are negative and the node is stable, thus Pis asymptotically stable. □ Electronic Research Archive Volume 32, Issue 2, 762–778. 768 2.2. Dynamics at the infinity In order to study the global dynamics of the system, we will use the Poincaré compactification, as it allows us to control the dynamics of a polynomial differential system near infinity. The following is an introduction to the basic notions of this technique; but more details can be found in Chapter 5 of [27]. Consider the sphere S2=ny=(y1,y2,y3)∈R3:y2 1+y2 2+y2 3=1o, which is called the Poincaré sphere. In general, a planar polynomial system of the form ˙x1=P(x1,x2), ˙x2=Q(x1,x2), can be projected into the sphere, thus obtaining an induced vector field in S2\S1. Note that we can identify R2with the tangent plane to the sphere at the point (0,0,1). Then, we can obtain the induced vector field by means of the central projections f+:R2→S2and f−:R2→S2, as defined by thw following: f+(x)= x1 ∆(x),x2 ∆(x),1 ∆(x)!and f−(x)= −x1 ∆(x),−x2 ∆(x),−1 ∆(x)!, where ∆(x)=qx2 1+x2 2+1. The differential D f +and D f −provide a vector field in the northern and southern hemispheres respectively, and we can analytically extend this vector field to the points of the equator by multiplying the field by yd 3, where dis the degree of the original vector field in R2. This is important as the points of the equator, S1of S2, correspond with the points of R2at infinity. To make calculations, we work in the local charts, (Ui, ϕi) and (Vi, ψi), of the sphere S2, where Ui=ny∈S2:yi>0o,Vi=ny∈S2:yi<0o,ϕi:Ui−→ R2, and ψi:Vi−→ R2for i=1,2,3, where ϕi(y)=ψi(y)=(ym/yi,yn/yi)for m<nand m,n,i. The extended field is called the Poincaré compactification of the original vector field. The expression of the Poincaré compactification in the local chart (U1, ϕ1) is as follows ˙u=vd"−u P 1 v,u v!+Q 1 v,u v!#,˙v=−vd+1P 1 v,u v!; (2.3) in the local chart (U2, ϕ2), the expression is as follows ˙u=vd"P u v,1 v!−uQ u v,1 v!#,˙v=−vd+1Q u v,1 v!,(2.4) and in the local chart (U3, ϕ3), the expression is as follows: ˙u=P(u,v),˙v=Q(u,v).(2.5) In the charts (Vi, ψi), where i=1,2,3, the system is the same as in the charts (Ui, ϕi) multiplied by (−1)d−1. As we want to study the behavior near infinity, we must study the infinite singular points (i.e., those which lie on the equator of the sphere). Electronic Research Archive Volume 32, Issue 2, 762–778. 769 It is sufficient to study the infinite points on the local chart U1and the origin of the local chart U2, because if y∈S1is an infinite singular point, then −yis also an infinite singular point; moreover, they have either the same or opposite stability depending on whether the system has odd or even degree. In the case of this work, since our system is motivated by a real population problem, it will only be of interest to study it for positive variables. Consequently, we shall present the phase portraits of the polynomial differential system in the positive quadrant of the Poincaré disc, which is the orthogonal projection of the northern hemisphere of S2 onto the plane y3=0. According to (2.3), in the chart U1, system (2.1) has the following expression: ˙u=−c1uv2+au2−(n+r)uv +c2v2+bu, ˙v=−c1v3+auv −rv2.(2.6) The singular points over v=0 are the origin of U1and the point (−b/a,0). As this second point is not on the positive quadrant of the Poincaré disk, it is not relevant in the problem we are studying. The linear part of system (2.6) at the origin is the following matrix: L= b0 0 0!, thus we have a semi-hyperbolic singular point. In order to determine its phase portrait, Theorem 2.19 in [27] can be used. As the singular point is the origin of the system and the Jacobian matrix is in the real normal form, we just compute the functions f(v) and g(v) in the theorem, thus obtaining the following: g(v)=−c2 bv2+o(x2). According to the theorem, it turns out that m=2 and am=−c2/b,0, and then the singular point is a saddle-node. Thus, we know that the sectorial decomposition of the singular point has two hyperbolic sectors and one parabolic sector, though it is necessary to determine the orientation and position of the three sectors. Studying the flow over the horizontal axes, we have that ˙v|v=0=0 and ˙u|v=0=au2+bu. Then, in a neighborhood of the origin, there are four possibilities for the position and orientation of the different sectors, which are given in Figure 2. u v (a) u v (b) v u (c) u v (d) Figure 2. Different configurations for a saddle-node in the origin of system (2.6) and (2.7). Electronic Research Archive Volume 32, Issue 2, 762–778. 770 Studying the flow over the vertical axis, we have that ˙u|u=0=c2v2>0 and ˙v|u=0=−v2(c1v+r). Note that this is not possible in the configurations (a)–(c) in Figure 2; therefore, the only possible phase portrait for the origin of the chart U1is the one given in Figure 2(d). According to (2.4), in chart U2, system (2.1) writes the following: ˙u=−c2uv2−bu2+(n+r)uv +c1v2−au, ˙v=−c2v3−buv +nv2.(2.7) The origin of this chart is a singular point, which is the only point at infinity that we have to consider here. The Jacobian matrix of the system at the origin is as follows: J= −a0 0 0!, where there is a semi-hyperbolic singular point. In this case, by applying Theorem 2.19 in [27], the following function is obtained: g(v)=c1/av2+o(x2); thus, according to the theorem, we have a saddle-node as m=2 and am=c1/a,0. Again, it is necessary to determine the position and orientation of the different sectors of the saddlenode. Analyzing the flow over the horizontal axis, we have that ˙v|v=0=0 and ˙u|v=0=−bu2−au. Then, in a neighborhood of the origin, there are four possibilities for the position and orientation of the different sectors, which are the same as those given in Figure 2. However, by studying the flow over the vertical axis, the only realizable configuration is the one in Figure 2(a), as we have ˙u|u=0=c1v2>0 and ˙v|u=0=v2(−c2v+n). 2.3. A note on the existence of a limit cycle The analytical proof of the existence of the limit cycle has turned out to be complicated. We have tried to apply Theorem 3.3. in [28], but although the bifurcation threshold with respect to certain parameters is easy to obtain, and the first genericity condition has been proven in some cases, tedious computations have led to a first Lyapunov coefficient equal to zero (see Appendix). The computations for the second Lyapunov coefficient were not manageable. Then, we carried out some numerical simulations that point to the fact that the limit cycle could appear under certain conditions. In Figure 3, we represent the solutions of system (1.1), with parameters r=a=b=n=0.2, c1=0.3, and c2=0.18. In this case, from an ecological point of view, after an increase and decrease in the number of individuals, the two species stabilize in the population corresponding to the equilibrium point; then, the number of prey and predators will be constant in the future time. In other words, the solutions tend to the equilibrium point of the coexistence. In Figure 4, the values of the parameters are r=50, a=2553/524288, c1=5, b=5/2147483, n=1000/1024, and c2=1. We observe an oscillatory behavior that could represent that of an orbit approaching the limit cycle. From an ecological point of view, this situation represents an oscillatory behavior where the number of prey and predators increase and decrease periodically. Electronic Research Archive Volume 32, Issue 2, 762–778. 777 as for the example q=(q1,q2)Tand p=(p1,p2)T, where q1=1 8bc2r (n−r)p4bc1+n2+n q2n−r2p4bc1+n2+n2p4bc1+n2+n3+nr2+4bc1(n+r)2+2n2!+ +2bc1 np4bc1+n2−3r−rp4bc1+n2+q2n−r2p4bc1+n2+n2p4bc1+n2+n3+nr2+4bc1(n+r)2+3n2+2r2!!, q2= np4bc1+n2+r−rp4bc1+n2−q2n−r2p4bc1+n2+n2p4bc1+n2+n3+nr2+4bc1(n+r)2+n2 4n−iω, p1=1 D bc2np4bc1+n2−3r−rp4bc1+n2−q2n−r2p4bc1+n2+n2p4bc1+n2+n3+nr2+4bc1(n+r)2+n2 2n(n−r)p4bc1+n2+n+2bc1, p2=− np4bc1+n2+r−rp4bc1+n2−q2n−r2p4bc1+n2+n2p4bc1+n2+n3+nr2+4bc1(n+r)2+n2 4nD −iω D, where ω=ω(a0) and D=1 4 rp4bc1+n2−np4bc1+n2−q2n−r2p4bc1+n2+n2p4bc1+n2+n3+nr2+4bc1(n+r)2−4bc1−n(n+r)!+ω2. Now, we compute the following: g20 =⟨p,B(q,q)⟩=g1 20 g2 20 , g11 =⟨p,B(q,q)⟩=−g1 11 g2 11 , g21 =⟨p,C(q,q,q)⟩=0, Electronic Research Archive Volume 32, Issue 2, 762–778. 778 where g1 20 =nr(2bc1+(n−r)(n+T)) 2bc13n2+n(T−3r)+2r2−rT +S+(n−r)2n2+S(n+T)  8r−2iv t16n2r2q4bc1nr +1 4n2+n(r+T)−rT2−1 4(n−r)2n2+T(n+r)−nr +S2 n2r2  −2in v t16n2r2q4bc1nr +1 4n2+n(r+T)−rT2−1 4(n−r)2n2+T(n+r)−nr +S2 n2r2+2n2+2n(r+T)−2(rT +S) , g2 20 =16c24b2c2 1n4+6n2r2+r4+n2(n−r)(n+T)2n2n2−2nr +5r2+S(n+r)(n−3r)+bc18n6+n5(−18r+4T) +nr2−Tq4bc1(n+r)2+2nn3+n2T+nr2−r2T+2r3+3rS + +3n416r2−2rT +S +r3Tq4bc1(n+r)2+2nn3+n2T+nr2−r2T−2r2T−2rS +n2r−Tq4bc1(n+r)2+2nn3+n2T+nr2−r2T +8r3−16r2T−5rS +n3Tq4bc1(n+r)2+2nn3+n2T+nr2−r2T−32r3+20r2T−7rS , g1 11 =abc1nr24n2+n(T−3r)−rT −S2b2c2 1+2bc1(n−r)(2n−r+T)+n(n−r)2(n+T) +(2bc1+(n−r)(n+T))2−2n2−2n(r+T)+2(rT +S) +i2nv t16n2r2q4bc1nr +1 4n2+n(r+T)−rT2−1 4(n−r)2n2+T(n+r)−nr +S2 n2r2  , g2 11 =c2nr n2+n(T−3r)−rT −S2bc13n2+n(T−3r)+2r2−rT +S+(n−r)2n2+S(n+T) +(2bc1+(n−r)(n+T)) 16n2r2r4bc1nr +1 4n2+n(r+T)−rT2−1 4(n−r)2n2+T(n+r)−nr +S2, and T=p4bc1+n2, S=q4bc1(n+r)2+2nn3+n2T+r2(n−T). Then, the first Lyapynov coefficient is as follows: ℓ1=1 2ω2Re(ig20g11 +ωg21)=0. ©2024 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0) Electronic Research Archive Volume 32, Issue 2, 762–778.