A prey-predator model with migrations and delays

In this paper we propose a prey-predator model in multiple 
patches through the stage structured maturation time delay with migrations among patches. Focus on the case with two patches, we discuss the existence of equilibrium points and the uniform persistence. In particular, when the maturation times are the same in the patches, we study the local and global attractivity of boundary equilibrium point with general migration function and the local stability of the positive equilibrium with constant migration rate. Numerical simulations are provided to demonstrate the theoretical results, to illustrate the effect of the maturation time, the migration rate on the dynamical behavior of the system.

(Communicated by Yang Kuang) Abstract. In this paper we propose a prey-predator model in multiple patches through the stage structured maturation time delay with migrations among patches. Focus on the case with two patches, we discuss the existence of equilibrium points and the uniform persistence. In particular, when the maturation times are the same in the patches, we study the local and global attractivity of boundary equilibrium point with general migration function and the local stability of the positive equilibrium with constant migration rate. Numerical simulations are provided to demonstrate the theoretical results, to illustrate the effect of the maturation time, the migration rate on the dynamical behavior of the system.

1.
Introduction. The prey-predator relationship is prevalent in the nature and hence it is one of the major themes in mathematical models. In the 1920s A. Lotka and V. Volterra formulated the original model to describe the prey-predator interactions. Since then, different prey-predator models have been built by incorporating additional biological processes into the classical LotkaVolterra model [2, 6, 8-10, 15, 17, 18, 20, 22, 23, 32, 36, 37]. For instance, in [8], the authors studied the dynamics of a diffusive prey-predator system with Allee effects in the prey growth. In [15], a prey-predator model with a consideration of the prey refuge has been studied. It is well-known that some kind of time delays are unavoidable in prey-predator population and they can affect the stability of equilibrium points [19]. In [6], the authors presented a stage-structured predator prey model with gestation delay. A delay prey-predator system without instantaneous negative feedback has been studied in [32]. The authors in [37] considered a class of delayed Lokta-Volterra prey-predator model with two delays. Time delay due to maturation of the predator is inescapable since the juvenile predators usually live on their parents until they become mature, or on a resource that is not wanted by adults and available in abundance. The effect of maturation delay on the behavior of populations has been studied in [12].
In prey-predator system with multi-patch environment, migration occurs when a species moves from its patch to another due to different factors like competition, age, sex, lack of food, climate and the season which make it one of the prevalent phenomena in nature [26]. For example, in aquatic environments, many zooplankton species exhibit vertical movements each day due to light and food. During the day time, some species migrate downwards into the darkness to reduce the risk of predation by fish, while at night time, these species move upward to consume the phytoplankton [4,24]. Empirical studies have shown that the interaction between patches has important consequences for population stability and persistence e.g. [11,31,33]. The influence of such interaction, migration, capture scientists' attention in the communities. Some researchers considered different migration rates for both predator and prey species. For instance, the authors constructed a two patches model with migration of both the predator and the prey in [10], and assumed the migration of the predator depend on the population of the prey in each patch, while the prey migration is in a constant rate. They found upper and lower bounds for the populations and discussed the stability and instability of the positive and boundary equilibrium points. In [17], a model was provided with two patches where the prey does not migrate and the predator in the patch with a higher density will migrate to the patch with a lower density. In [22], it was assumed that the predator migrates with a constant migration rate, while the prey migration depends on the predator density, the authors proved the existence of a unique strictly positive equilibrium. In [2], it was assumed that the prey (predator) migration rate depend on the predator (prey) density and showed that for a large class of density-dependent migration rules there exists a unique and stable equilibrium for migration.
To the best of our knowledge, there is not so many result in the literature in the research of joined consideration of migration and time delay in biological models. Introducing time delays to prey-predator models usually brings challenges to the mathematical analysis and may change the dynamics. In [36], a model with two prey and one predator population, time delay due to gestation and constant prey migration rate was investigated. The authors discussed the persistence, derived conditions for the existence and local/global stability of the positive equilibrium. All the models studied in [2,10,17,22] have not considered the influence of time delay in any evolution stage although migrations were collaborated. In this paper, we propose a general model with multiple parallel patches which combined the migration interactions among all patches and constant maturation time delay through stage structure. Following the biological feasibility, we take general functional responses with respect to migration rates in both prey and predator classes to allow a unified treatment for all biological important cases. From the viewpoint of dynamics, we analyze the model with proving the positivity and boundedness of the solutions, local stability of the steady states, global attractivity of boundary equilibrium point and uniform persistence of the system. Then, we give numerical simulations to show the theoretical results and to see how the dynamical behavior is sensitive to the maturation time of predator, the migration rates in prey and predator.
The paper is organized as follows: in Section 2, a general model of multiple parallel prey-predator patches with migration of both species is studied. When the model contains two patches, we provide conditions for the uniform persistence, discuss the existence of boundary equilibrium points, and further investigate their stabilities with the same time delays. In Section 3, we explore the existence and stability of positive equilibrium point for a special case when the two patches are symmetric and the prey and predator migrations are constants. Numerical simulations are given to compliment the theoretical analysis for such a general system in Section 4. Conclusion is drawn in Section 5.

2.
A model of multiple parallel prey-predator patches with migration and maturation delays. In previous work [3], we studied a model including multiple parallel food chains, each consisting of a prey species P i and its dedicated predator Z i , with maturation delay and without interaction between patches. In this paper we consider the interaction between patches, assume all P i compete for limiting nutrient, only adult predators are capable of preying on the prey species and propose a more complex model for the multi-patch prey-predator system, with prey (predator) migration rates depending on the corresponding predator (prey) density: where i = 1, . . . , n, µ i (N ) is the growth rate of P i as a function of nutrient concentration N = T − n k=1 α k P k − n k=1 β k Z k ≥ 0, with the coefficients α k , β k , (k = 1, · · · , n) related to the efficiency of nutrient consuming for each species. The author in [5] used the form N = T − n i=1 P i − n i=1 Z i which treats the consuming abilities for all the species equally. As we know that in natural world, different species may have different strength in taking nutrient, so to be more realistic, we give different weight on nutrient consumption for the species. The constant parameter T is the internal nutrient supply, h i (P i ) is the per-prey-per-predator harvest rate of P i by Z i and i n j=1 Z j is predation on the ith predator by higher trophic levels which is assumed to be proportional to the total predator n j=1 Z j (see [5]); b i denotes the adult predator's birth rate, d i is the mortality death rate of the juvenile (through-stage death rate), τ i is time to mature and f ij (Z i ) (g ij (P i )), j = 1, . . . , n and i = j, represents the migration rate for prey (predator) from patch i to patch j. All the parameters are positive constants. An architecture of the model (1) with 3 patches is given in Fig. 1.
From the biological point of view, we assume all the functions in (1) are continuous and differentiable over R + , and satisfy the following hypotheses for i, j = 1, . . . , n with i = j: . The assumptions of (C 1 ), (C 2 ) and (C 3 ) are biologically obvious, (C 4 ) and (C 5 ) indicate that the preys escape from a patch more actively when more predators are in the same patch, while if more preys stay in a patch, the movement of the predator to other patch is inactive .
The general model (1) can cover several general, partial general and specific models in literature with some particular choice of the functions. To name only a few models with two patches without time delay, the system in [20] is a special model of (1) with µ(N ) = µ(P ), h(P ) = p(P ) P , same , f (Z) = constant and g(P ) = 0; when we take µ(N ) = r(1 − P K ), h(P ) = b b+P , g(P ) = d 1+P , (Z) and f (Z) are constant, the model (1) is the same as that in [10]; When one choose f ij (Z) = β i Z + β 0 and g ij as constant, the migration terms in (1) become the fast model in [22]; while we restrict f and g with df dZ > 0 and dg dP < 0, then the migration terms in (1) become the fast model in [2].
Let τ = max{τ i , i = 1, . . . , n} and define C : Then C is a Banach space and C + = C([−τ, 0), R 2n + ) is a normal cone of C with nonempty interior in C. We have the following result for the positivity and boundedness properties of the system (1).
With general n, the idea is similar by using mathematical induction.
To prove the uniqueness of suchẼ, it's easy to see that when one of P i , (i = 1, 2) is zero, another one must be zero as well, which excludes the case with one prey species can survive and another cannot.

Remark 2.
With migrations between patches, only one positive predator-free equilibrium point exists for system (1). But without migration, the system has infinite number of predator-free equilibrium points with the form [3]. The result in Proposition 1 implies the effect of migrations.
To study the stability of each equilibrium point (x 1 ,y 1 ,x 2 ,y 2 , . . . ,x n ,y n ), we can obtain the general form of the characteristic equation, which is, Due to the complexity of the construction in the characteristic equation with delay, we only discuss the stability of boundary equilibrium points when n = 2.
Theorem 2.2. E 0 is always an unstable steady state.
Proof. The linearization of (2) at E 0 is The characteristic equation of (4) is which is equivalent to Thus 2µ1(T )µ2(T ) (µ1(T )+µ2(T )) 2 > 1 leads to a contradiction. Therefore h(λ) = 0 has at least one positive real root, implying E 0 is unstable steady state.
It is easy to show that all the roots of det(B) = 0 are negative real roots by using Gershgorin circle theorem in [27] and the symmetry of To discuss the stability of the predator-free equilibrium pointẼ = (T , 0, T , 0), we introduce the following result given in [19] to discuss the global attractivity of E. Lemma 2.3. In a delay system where a, b, τ > 0, the following results hold For the special case with τ 1 = τ 2 := τ , we knowẼ is globally attractive under certain condition, which is given in the following theorem.
In the stability analysis, we are more interested in the influence of the key parameters, functions including the migration rates f ij , g ij and the time delays τ i (i = 1, 2). We notice that K(λ, τ 1 , τ 2 ) is independent of the prey migration rates f ij , so the prey migrations do not affect the stability of the predator-free equilibrium pointẼ. On the other hand, if there is no predator migration, i.e. g ij = 0, [3]; but with general predator migration, the stability atẼ becomes complicated. In the following, we focus on the stability with respect to time delay only and we will investigate the influence of migration rates numerically.
Under the assumption thatẼ is locally asymptotically stable at τ * , we have that E is locally asymptotically stable on this interval, since there are no sign changes in roots of K(λ, τ ) = 0 for τ ∈ (τ,τ ).
Adapting ideas from [35,38], we can study the persistence of the system (2), when one of the following assumptions holds: The following result shows that system (2) is uniformly persistent if either (R 1 ) or (R 2 ) holds.
Proof. We prove the theorem under the first assumption of (R 1 ). For the second case, a similar proof works. Define and Claim 1. There exists a δ 1 > 0, such that for any φ ∈ X 0 , By contradiction, suppose that lim sup By (C 4 ), if φ 1 (0) > 0 and φ 3 (0) > 0, then
Without loss of generality, we assume that (R 1 ) holds. Choose ζ > 0 small enough such that Let t ζ > t 1 + τ be sufficiently large such that Thus, by (C 3 ) and (C 5 ) for t > t ζ By the comparison principle and Lemma 2.3, we get Z 1 (t) → ∞ as t → ∞, which is a contradiction.
Define a continuous function p : X → R + by It is clear that p −1 (0, ∞) ⊂ X 0 and if p(φ) > 0 then p(Φ(t)φ) > 0 for all t > 0. By Claim 3, we get that for any forward orbit of Φ(t) in M ∂ converges to E 0 orẼ, by Claim 1 and 2, we conclude that E 0 andẼ are two isolated invariant in X, and (W s (E 0 ) ∪ W s (Ẽ)) ∩ X 0 = ∅ and no subset of {E 0 ,Ẽ} form a cycle in ∂X. By [30], it then follows that there exists η > 0 such that lim inf t→∞ p(Φ(t)φ) ≥ η for all φ ∈ X 0 , which implies the uniform persistence. This completes the proof. 3. With constant migration rates. After we obtain the conditions for the existence and stability of the boundary equilibrium points E 0 andẼ, we simplify the model (1) and assume that all the prey (predictor) migration rates are the same constant f ij (Z i ) = m (g ij (P i ) = s) to explore the existence and local stability of positive equilibrium solutions for i = 1, 2, j = 3 − i. With this simplification, the system (2) becomes When the system (16) exists a positive equilibrium point E * = (P * 1 , Z * 1 , P * 2 , Z * 2 ), then P * i , Z * i , i = 1, 2 must satisfy with Otherwise, it leads to a contradiction with the summation of the first and third equations in (17). From the second and forth equations of (17), we have .
The condition (H 2 ) is necessary for the existence of a positive equilibrium point E * = (P * 1 , Z * 1 , P * 2 , Z * 2 ) in (17), but is not a sufficient condition in general. Define H i (P i ) = P i h i (P i ) for i = 1, 2. From the hypothesis (C 2 ), H i (P i ) is positive and increases in (0, T α1+α2 ). We can obtain from the second and forth equations of (17), respectively, where Z * 1 and Z * 2 are the implicit positive solutions of the following equations: which is impossible to find an analytical solution in general.
With fixed parameters, feasible values of delays and certain functions, for instance, choosing the functions given in Table ( When the two parallel food chains are symmetric, that is, the growth for all patches have the same functional forms, parameter values and the maturation delay, the system (16) becomes . Let τ s = sup {τ ∈ (0, τ max ) |system (19) has at least a positive equilibrium point } .
When the per-prey-per-predator harvest rate h(P ) is increasing, we know the positive equilibrium point E * exists with the form of (P * , Z * , P * , Z * ).
Proposition 3. Assume h(P ) is an increasing function on the interval (0, T 2α ). If the positive equilibrium point E * exists, then it must have the form of (P * , Z * , P * , Z * ).
In the following example, we consider a symmetric system with constant migration rates. First, we fix the value of the migration rates and study the dynamic behavior of system (19) as the time delay varies. Then, we fix the time delay and one of the migration rates to explore the affection of another migration rate on the stability of the positive steady state.
At τ = 0, E * (0) = (1.22,1.47,1.22,1.47) is stable by checking the condition (H ΦΨ 0 ). With τ >0 and fixed values g 12 (P 1 ) = g 21 (P 2 ) = 0.1 and f 12 (Z 1 ) = f 21 (Z 2 ) = 0.6, we plot S Ψ 0 (τ ), S Ψ 1 (τ ) and S Ψ 2 (τ ) defined in Eq. (24) in Fig. 2. We can see that the curve of S Ψ 1 (τ ) is negative when τ < τ * 1 = 0.41 and then becomes positive as τ increases up to τ * 2 = 24.11; and return to negative after that, implying for small delay τ < τ * 1 , the positive equilibrium point E * (τ ) is stable and loses the stability when τ ∈ (τ * 1 , τ * 2 ) and back to stable for large time delay, i.e τ > τ * 2 , before it disappears at τ = 45.02. These tell us the occurrence of stability switch, which is confirmed in Fig. 3. Biologically, when the maturation process is too short, both of the prey and the predator move to a certain level with constant concentration; while the maturation time increases, all the species oscillate regularly. Nevertheless, when the predator is mature enough, it behaviors less active which gives a chance to prey to increase. Therefore, the state of coexistence of the predator and the prey is stable, rather than oscillatory.
We notice that the chosen function h(P ) in (31) is not a increasing function, from Remark 4 it is possible to have other positive equilibrium point beside E * (τ ) but we failed to find it although we have tried different choices of the functions and parameters.
Next, we discuss how the variation of prey and predator migration rates effect the dynamical system. Since the system is symmetric, let f 12 (Z 1 ) = f 21 (Z 2 ) = m, g 12 (P 1 ) = g 21 (P 2 ) = s. In general, the stability is corresponded to the roots in Φ(λ, τ ) = 0 given in (25), where Φ depends on both migration rates m, s, and Ψ(λ, τ ) = 0 given in (26)   In the second example, we consider a non-symmetric system with constant and non-constant migration rates.
Example 2. A non-symmetric system. Choose First, we take a constant migration rates f 12 (Z 1 ) = f 21 (Z 2 ) = m and g 12 (P 1 ) = g 21 (P 2 ) = s. With fixed m = 0.6 and s = 0.8, in Fig. 5, we can observe interesting dynamical behavior as τ varies. When the maturation time is τ = 5, a stable steady state exists (Fig. 5a); with the increasing of the maturation time, the steady state loses the stability and the system becomes oscillatory (Fig. 5b); With further increasing of τ , it is interested to observe the occurrence of period-doubling (Fig.  5c), and then regain the stability (Fig. 5d) when τ is even larger.
To see the influence of the prey migration rate m, we fix s = 0.8 and τ = 5. Fig.  6 shows that the increasing of m can destabilize the steady state. Corresponding to the remarkable behavior in Fig. 5c, where m = 0.6, s = 0.8, τ = 20, we keep all the parameters and delay the same, just slightly change the value of m around m = 0.6,            Similarly, we fix m = 0.6 and vary the value of s. When s = 0.6, there exists periodic oscillation. Increasing s from 0.6 to 0.7 or 0.8, another oscillatory behavior appears to form a periodic-doubling. It is very interesting to observe that the stability switch occurs with the further increasing of s, see Fig. 8.
If we choose the migration rates in (31) as non-constant functions, where ϑ is a positive real number. First, we fix ϑ = 0.5. Fig. 9 shows the occurrence of stability switch and a periodic-doubling when the time delay varies. Then we pick the moment shown in Fig. 9c, change the value of ϑ slightly around ϑ = 0.5. It is surprising that the periodic-doubling is broken when we increase or decrease ϑ, see Fig. 10, implying that the system may have rich dynamics with the change of migration rate.
In the last example, we consider a non-symmetric system having three patches and with constant migration rates. We compare this system with the model without migration rates which is studied in [3].     . Time series P 2 (t) and Z 2 (t) of system (2) with parameters and functions given in (30), (32) and (33).
If there is no migration, i.e. m = s = 0, when the maturation time is very small, a doubly periodic solution exists [3], see Figs. 11. While when m and s are very small, then the doubly periodic orbit is disappeared, although a periodic solution still exists, see Fig. 12a. However, when either m or s is relative large, the periodic solution is vanished and the system becomes stable, see Fig. 12b-12d. Therefore the dynamical behavior is sensitive to prey and predator migration rates. Seems the active migration for prey or predator can stabilize the ecosystem. 5. Conclusion. In this paper, we propose a general model with multiple parallel food chains through the stage structured maturation time delay, and consider      Figure 11. Phase portrait of (16) with parameters and functions given in (30) and (34), m = s = 0.  Figure 12. Phase portrait of (16) with parameters and functions given in (30) and (34).
predator and prey migrations among all patches. We have carried out mathematical analysis to discuss the existence of the steady states and their stabilities. Taking the system with two patches as an example, we have proved the existence and uniqueness of two boundary equilibrium point E 0 and the predator-free equilibrium pointẼ. We have shown that E 0 is always unstable, provided conditions for the uniform persistence, obtained sufficient condition for the global attractivity ofẼ and discussed its local stability with the same time delays. When both patches are symmetric with constant migration rates, we have discussed the existences and local stability of positive equilibrium point; obtained that the time delay can not only destroy the existence, but also destabilize the positive equilibrium even it exists. To complement the analytical results, we have used numerical techniques to examine how the maturation time of predator, the prey and predator migrations affect on the dynamical behavior. The numerical simulations show the rich dynamics including stability switch and periodic-doubling.