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

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,4,5,6,7,8,9,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 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,15,16,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 Ⅰ 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 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:

where $r$ represents the growth rate of prey, $n$ is the death rate of predator, $a$ is the rate of predation, $b$ is the conversion rate of eaten prey into new predators, $h$ and $\alpha$ are the functional response coefficients (which involve for example the handling time), and $c_1$ and $c_2$ are 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 $\alpha = h = 0$, which corresponds with a Holling type Ⅰ functional response. Our main result is as follows.

Theorem 1.1. The global phase portrait of system (1.1) where $\alpha = h = 0$, in the positive quadrant of the Poincaré disc, is one of the following:

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 = \alpha = 0$, corresponds to a Holling type Ⅰ functional response:

In order to make the reading more clear, we introduce the following notation:

There are two general solutions for $\dot{x} = \dot{y} = 0$, with the following expressions

which are always well defined since the parameters are positive and $R$ is 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 = \left\lbrace(x, y) \in \mathbb{R}^2 \mid x \geq 0, y \geq 0 \right\rbrace$ is positively invariant.

Proof. The fact that region $C$ is positively invariant can be deduced from the fact that

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 $\left\lbrace(0, y) \in \mathbb{R}^2 \mid y \geq 0 \right\rbrace$ and $\left\lbrace(x, 0) \in \mathbb{R}^2 \mid x \geq 0 \right\rbrace$ are 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

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, $-b c_1 - a c_2 + nr + R < 0$. Then,

by squaring the expressions,

which contradicts the fact that all parameters are positive. An analogous reasoning allows us to prove that the other solution of $\dot{x} = \dot{y} = 0$ is never positive. If we suppose that it is positive, then $-b c_1 - a c_2 + nr - R > 0$; by squaring the expressions, one can obtain the following:

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 $\left((n+r)(b c_1 + a c_2) + (n-r)(nr + R)\right)^2 < 16 n^2 r^2 R$.

2) It is a stable node if $\left((n+r)(b c_1 + a c_2) + (n-r)(nr + R)\right)^2 > 16 n^2 r^2 R$.

Proof. The Jacobian matrix of system (2.1) at the singular point $P$ is

and has the following eigenvalues:

Then, the trace of the Jacobian matrix is as follows:

where the expression is always negative, and can be proved. First,

is negative as

Furthermore, taking into account that

it holds that

which is negative as $b c_1 + a c_2 + n r > \sqrt{(b c_1 + a c_2 + nr)^2-4 n r a c_2}.$ Then, we have that the trace is always negative for any values of the parameters.

The determinant of the Jacobian matrix is as follows:

which is positive.

If the radicand in the expressions of $\lambda_1, \lambda_2$ is negative (i.e., if $\left((n+r)(b c_1 + a c_2) + (n-r)(nr + R)\right)^2 < 16 n^2 r^2 R$), then the two eigenvalues are complex, and the singular point $P$ is a stable focus as $\lambda_1+\lambda_2 = 2 \text{Re}(\lambda_1) = 2 \text{Re}(\lambda_2) < 0$. Then, $P$ is asymptotically stable.

If the radicand is positive (i.e., if $\left((n+r)(b c_1 + a c_2) + (n-r)(nr + R)\right)^2 > 16 n^2 r^2 R$), then the eigenvalues are real; and since $\lambda_1\lambda_2 > 0$, the singular point is a node, and taking into account that $\lambda_1+\lambda_2 < 0$, then both eigenvalues are negative and the node is stable, thus $P$ is asymptotically stable.

□

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 $\mathbb{S}^2 = \left\lbrace y = (y_1, y_2, y_3) \in \mathbb{R}^3 : y_1^2+y_2^2+y_3^2 = 1 \right\rbrace$, which is called the Poincaré sphere. In general, a planar polynomial system of the form

can be projected into the sphere, thus obtaining an induced vector field in $\mathbb{S}^2 \backslash \mathbb{S}^1$. Note that we can identify $\mathbb{R}^2$ with 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^+:\mathbb{R}^2\to \mathbb{S}^2$ and $f^-:\mathbb{R}^2\to \mathbb{S}^2$, as defined by thw following:

where $\Delta(x) = \sqrt{x_1^2+x_2^2+1}$.

The differential $Df^+$ and $Df^-$ 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 $y_3^d$, where $d$ is the degree of the original vector field in $\mathbb{R}^2$. This is important as the points of the equator, $\mathbb{S}^1$ of $\mathbb{S}^2$, correspond with the points of $\mathbb{R}^2$ at infinity.

To make calculations, we work in the local charts, $(U_i, \phi_i)$ and $(V_i, \psi_i)$, of the sphere $\mathbb{S}^2$, where $U_i = \left\lbrace y \in \mathbb{S}^2 :y_i > 0 \right\rbrace$, $V_i = \left\lbrace y \in \mathbb{S}^2 : y_i < 0 \right\rbrace$, $\phi_i: U_i \longrightarrow \mathbb{R}^2$, and $\psi_i: V_i \longrightarrow \mathbb{R}^2$ for $i = 1, 2, 3$, where $\phi_i(y) = \psi_i(y) = \left(y_m/y_i, y_n/y_i \right)$ for $m < n$ and $m, n \neq i$.

The extended field is called the Poincaré compactification of the original vector field. The expression of the Poincaré compactification in the local chart $(U_1, \phi_1)$ is as follows

in the local chart $(U_2, \phi_2)$, the expression is as follows

and in the local chart $(U_3, \phi_3)$, the expression is as follows:

In the charts $(V_i, \psi_i)$, where $i = 1, 2, 3$, the system is the same as in the charts $(U_i, \phi_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).

It is sufficient to study the infinite points on the local chart $U_1$ and the origin of the local chart $U_2$, because if $y\in \mathbb{S}^1$ is an infinite singular point, then $-y$ is 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 $\mathbb{S}^2$ onto the plane $y_3 = 0$.

According to (2.3), in the chart $U_1$, system (2.1) has the following expression:

The singular points over $v = 0$ are the origin of $U_1$ and 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:

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:

According to the theorem, it turns out that $m = 2$ and $a_m = -c_2/b\neq0$, 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 $\dot{v}\mid_{v = 0} = 0$ and $\dot{u}\mid_{v = 0} = au^2+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.

Studying the flow over the vertical axis, we have that $\dot{u}\mid_{u = 0} = c_2v^2 > 0$ and $\dot{v}\mid_{u = 0} = -v^2(c_1v+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 $U_1$ is the one given in Figure 2(d).

According to (2.4), in chart $U_2$, system (2.1) writes the following:

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:

where there is a semi-hyperbolic singular point. In this case, by applying Theorem 2.19 in ^[27], the following function is obtained:

thus, according to the theorem, we have a saddle-node as $m = 2$ and $a_m = c_1/a\neq0$.

Again, it is necessary to determine the position and orientation of the different sectors of the saddle-node. Analyzing the flow over the horizontal axis, we have that $\dot{v}\mid_{v = 0} = 0$ and $\dot{u}\mid_{v = 0} = -bu^2-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 $\dot{u}\mid_{u = 0} = c_1v^2 > 0$ and $\dot{v}\mid_{u = 0} = v^2(-c_2v+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$, $c_1 = 0.3$, and $c_2 = 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$, $c_1 = 5$, $b = 5/2147483$, $n = 1000/1024$, and $c_2 = 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.

2.4.   Global phase portrait

To determine the global phase portrait on the positive quadrant of the Poincaré disk, we have to gather the local information obtained.

From Theorem 2.3 and 2.4, we know that there is always one positive singular point, and it is either a stable node or a stable focus, which are topologically equivalent.

At infinity, there are two singular points: the origins of the charts $U_1$ and $U_2$. The origin of $U_1$ has the phase portrait in Figure 2(d); as {the region of interest in our model is the} positive quadrant, there are orbits that enter into this quadrant from each point of the axis $u = 0$, except from the origin.

The origin of $U_2$ has the phase portrait in Figure 2(a). Again, focusing on the positive quadrant, a separatrix exists which leaves from this infinity singular point and gives rise to two different sectors, one between the separatrix and infinity and the other between the separatrix and the $u = 0$ axis. In the first one, there are orbits which all leave from the origin of $U_2$; in the second one, there are orbits that enter the positive quadrant of the Poincaré disk from each one of the points of the $u = 0$ axis.

Accordingly, two global phase portraits can be obtained (i.e., those given in Figure 1), depending on whether or not there is a limit cycle surrounding the positive singular point.

In case there is such a limit cycle, the separatrix, as well as all orbits entering the positive quadrant from each of the points on the $u = 0$ and $v = 0$ axes, go to the limit cycle when $t$ tends to infinity. In case it exists, the limit cycle is an attractor on the outside and a repulsor on the inside, and all orbits starting from some point in the region enclosed by the limit cycle go to the positive singular point when $t$ tends to infinity.

If there is no limit cycle, all orbits go to the positive singular point, which will then be globally asymptotically stable (in the positive quadrant of the Poincaré disk).

3.   Discussion and conclusions

In this work, we have studied a predator-prey system with a Holling type Ⅰ functional response, with the particularity that we have added constant terms representing immigration in both the predator and the prey species.

The main result of our work is that two dynamical behaviors can appear in a predator-prey ecosystem with this characteristic, which are given in Figure 1.

We highlight an important outcome: in both cases, there is an asymptotically stable singular point, which, from a biological point of view, means that both populations can coexist. In the first case, the singular point is globally asymptotically stable, thus the result is even stronger: regardless of the initial number of prey and predators, both species tend to coexist. If a limit cycle exists, then it delimits the basins of attraction of the asymptotically stable singular point; then, if the initial condition of the number of prey and predators is inside the region delimited by the limit cycle, then the populations also tend to coexist. In the other case, for the rest of the initial conditions, the trajectories tend to the limit cycle, which is an attractor on the outside; thus the populations tend to a cyclical behavior. These results are topologically represented in Figure 1.

It has not been possible to analytically prove the existence of the limit cycle, which is why we have relied on some numerical simulations; for example, in Figure 4, these simulations allow us to see a certain oscillatory behavior, in which the number of prey and predators continually increases and decreases over time. Note that this does not allow us to definitively conclude the existence of the limit cycle.

These results differ from those obtained in the classical Lotka-Volterra model, in which there is a stable singular point that is topologically a center and all other solutions are periodic orbits. Thus, from our study, it can be concluded that the presence of a certain number of immigrants in both species, which is something that is not strange in the nature, can affect the dynamics of the ecosystem, thus leading to an asymptotically stable coexistence equilibria. The results point in the same direction as previous studies, in which immigration was considered only in one of the species, as in ^[20] or ^[22]. It is difficult to compare our results with some other works, where we do not see clearly the effect of immigration, as it appears to be combined with different types of properties of the ecological systems, as is the case of ^[21]. We have not found works in which immigration in both species is considered, nor works in which the global dynamics is studied; therefore we can compare with our results models that have conditions similar to ours.

As open problems for the future, it will be interesting to find a specific set of parameters for which the limit cycle exists, and also replace the constant immigration rate by other functions that can also represent the number of immigrants in both species.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The author is partially supported by the Ministerio de Ciencia e Innovación, Agencia Estatal de Investigación (Spain), Grant PID2020-115155GB-I00 and the Consellería de Educación, Universidade e Formación Profesional (Xunta de Galicia), Grant ED431C 2023/31 with FEDER funds.

Conflict of interest

The authors declare there is no conflicts of interest.

Appendix

As stated in subsection 2.3, by trying to prove the existence of a limit cycle that appears by a Hopf bifurcation, we have computed the first Lyapunov cofficient, although it turned out to be zero, so it was not possible to conclude the existence of the bifurcation. Here, we summarize the computations.

The Jacobian matrix at the equilibrium $P$ is as follows:

and has the eigenvalues $\mu(a)\pm \omega(a) i$, where

We obtain $\mu(a_0) = 0$ for the following:

Then, if the condition $16 n^2 r^2 R - \left((n+r)(b c_1 + a c_2) + (n-r)(nr - R)\right)^2 > 0$ holds, at $a = a_0$, the equilibrium point $P$ has a pair of pure imaginary eigenvalues $\pm i \omega (a)$ and the system will have a Hopf bifurcation if some Lyapunov constant is nonzero and $(d\mu/da)(a_0)\ne 0$. Note that we would be interested in this case when the expression for $a_0$ is positive, though this is not a problem.

It is necessary to check if the genericity conditions are satisfied (see ^[28,Theorem 3.3]). We check that the transversality condition is satisfied as

and with the help of the softwate Mathematica, v. 13.3.1, we have proven that this expression is always negative.

To check the second condition, the first Lyapunov constant must be calculated. We fix the value $a = a_0$; then, the equilibrium $P$ has the expression $P_2 = (P_{21}, P_{22})$, with

We translate $P$ to the origin of coordinates, thus obtainig a system that can be represented as follows:

where $A = A(a_0)$ and the multilinear functions $B$ and $C$, which are given by the following:

We need to find two eigenvectors $p, q$ of the matrix $A$ verifying

as for the example $q = (q_1, q_2)^T$ and $p = (p_1, p_2)^T$, where

where $\omega = \omega(a_0)$ and

Now, we compute the following:

where

and

Then, the first Lyapynov coefficient is as follows: