1 Introduction
The role of spatial heterogeneity and migrations on predator–prey or host–parasite systems stability has been studied in many papers. In [1], the authors exhibit the stabilizing effect of migrations in a spatially distributed host–parasite system. Some authors have then investigated the impact of density dependence in the migration processes on the predator–prey interactions stability (see [2–7], for instance). Even if the initial works have shown that the spatial heterogeneity and the migrations should stabilize the systems, Murdoch and Oaten [8] studied a system exhibiting the opposite effect. It followed a series of papers dealing with the impact of spatial heterogeneity on predator–prey (or host–parasitoid) systems' stability and the main conclusion was that the stability of the system depends on the details included in the models.
In the previous works, the interactions models had either asymptotically stable or unstable equilibria. For instance, in [7], the model is a two patches host–parasitoid system with a linear growth function for the host on each patch, a Holling type-II functional response on each patch and a linear decay for the parasitoid in the absence of host on each patch. In this kind of system, the positive equilibrium is known to be unconditionnally unstable. The density-dependent migration rates used by the author can stabilize the equilibrium. However, it is difficult to quantify exactly which part of the density dependence in the rates, which part of the spatial heterogeneity, and which part of the interaction processes contribute to the stabilization. In the present work, we expect that the choice of a neutrally stable model for the interaction part and constant migration rates should give a best insight into the understanding of the role of spatial heterogeneity on stabilization of predator–prey models.
Indeed, we consider a predator–prey model on two patches. One patch is a refuge for the prey, thus there is no predator in this patch. On the second patch, we used a Lotka–Volterra model to describe the interactions. This choice is governed by the stability property of this model: the positive equilibrium is neutrally stable. In this case, the interactions part do not contribute neither to stabilization nor to unstabilization. For the migrations, we used constant rates in order to eliminate any effect of density dependence. Since the interaction part is neutral from the stability point of view, all the stabilization effects are contained only in the migration processes and in the spatial heterogeneity. Moreover, we can easily admit that a good choice of the density dependence function used for the migration rates can compensate the terms generating the unstability. Since we have chosen density-independent migration rates, the stabilization is just the result of spatial heterogeneity.
In order to deal with the three-dimensional system, we use a time scale arguments which allows us to consider only the total populations densities, that is a two-dimensional system. This kind of methods has been described for instance in [9–11]. Our approach can be compared in a first approximation to the so-called ‘quasi-steady-state’ assumption. However, we show in the particular model of this paper that the later method is not sufficient and that the former can be used. The obtained two-dimensional system is a perturbation of the Lotka–Volterra model. This famous model is not structurally stable and our analysis shows that the perturbation ‘breaks’ the closed curves. In fact, for the complete initial model, there exists a two-dimensional invariant attracting manifold. The aggregated model is the restriction of the complete system to this invariant manifold. In practice, the invariant manifold is obtained by the use of a Taylor expansion with respect to a small parameter (the time scale factor). In the first approximation, the reduced model is of Lotka–Volterra type. Thus we need to continue the analysis by computing the second term in the expansion.
Finally, we prove that there is a unique equilibrium with positive coordinates and that this equilibrium is globally asymptotically stable, for all the parameter values. Thus the spatial heterogeneity has a stabilizing effects on the predator–prey interactions. Furthermore, there is no possible bifurcation on the invariant manifold, thus the previous result is complete.
The plan of this article is as follows: we first describe the complete predator–prey model, the following section explain how the reduced model is obtained. Then we analyse the reduced model and the last section provides a discussion on the results.
2 The complete model
We consider a predator–prey model on two patches. The prey can move on both patches while the predator remains on patch 1. The patch 2 is thus a refuge for the prey. We denote by the prey density on patch i, . We denote by y the predator density. On each patch, the prey population growth rate and the predator population death rate are linear, the predation rate is bilinear, that is proportional to prey and predator densities and the predator growth rate is proportional to the predation rate. The model is given by the following set of three ordinary differential equations:
(1a) |
(1b) |
(1c) |
Let be the total amount of prey population. It follows that and are the proportions of prey on patch 1 and patch 2, respectively. With these variables, we can write the system (1) in the following equivalent way:
(2a) |
(2b) |
(2c) |
3 Reduction of the dimension
In this section, we build a two-dimensional system governing the dynamics of the total populations densities x and y. Morevover, this system gives the same dynamics as that obtained for x and y in the system (2). This will facilitate the mathematical study of system (1) provided in the next section.
The reduction method is based on the fact that the system (2) has two different typical time scales. From the mathematical point of view, the method is based on a theorem due separately to Hirsch, Pugh and Shub in [12] and Fenichel in [13]. The interested reader can find more informations in [14,15] or in [16], for instance, and in the references given within these works. The mathematical theorems provide a justification to a method that is often used on the basis of intuition: since variables are fast, after a short time, these variables are close to their equilibrium values (if they exist) and then we replace the fast variables by their equilibrium values in the differential equations governing the slow variables, leading to a differential system that has the same dimension as the number of slow variables. This method is sometimes called ‘quasi-steady-state assumption’ method and can be used to build trophic chain models [17] and to analyse them [18,19]. However, in some cases, like in the present work, the quasi-steady-state assumption is not sufficient to determine the dynamics of the slow variables and the mathematical theorems provide usefull results that permit to conclude on the dynamics. We will illustrate this on the model studied in this paper.
Let us start to calculate the fast equilibrium, that is the equilibrium value of the fast variable . In order to get this equilibrium value, we put in system (2). The result is:
(3) |
(4a) |
(4b) |
The system (4) is a classical Lotka–Volterra model. All the solutions of this system with initial conditions in the positive quadrant are periodic. There is a positive equilibrium which is a centre. However, the dynamics of x and y in the system (2) does not match with the Lotka–Volterra dynamics, as illustrated in Fig. 1. Indeed, when we replace the fast variable by its equilibrium value, we make an error in the order of ɛ. Since the Lotka–Volterra model is not structurally stable, the ɛ-error can play an important role in the dynamics. The Fenichel theorem, for instance, claims that there is an invariant manifold defined by in the phase space . Since the fast equilibrium is hyperbolically stable, the manifold is normally hyperbolically stable, that is a trajectory starting in the neigbourhood of is attracted with an exponential speed. The previous approximation we made is thus a zero-order approximation of this manifold.
![](/biologies/article/CRBIOL_2004__327_11_1058_0/jpg/src/tex/figures/gr001.jpg)
Comparison between the dynamics of x and y given by the complete system (1) and that obtained with the two-dimensional system (4). We see on this figure that the complete system exhibits a stable focus while with the reduced system, the trajectory is a closed curve. The parameters values used in the simulation are: , , , , a=1, d=2, b=0.9 and ɛ=0.05.
We can get a first-order approximation of the manifold in the following way. Let us write:
(5) |
(6) |
(7) |
(8) |
(9a) |
(9b) |
![](/biologies/article/CRBIOL_2004__327_11_1058_0/jpg/src/tex/figures/gr002.jpg)
Comparison between the dynamics of x and y given by the complete system (1) and that obtained with the two-dimensional system (9). We see on this figure that the reduced system solution matches quite well the complete model solution. The parameters values used in the simulation are: , , , , a=1, d=2, b=0.9 and ɛ=0.05.
4 Mathematical analysis
In this section, the mathematical analysis of the system (9) is performed. We show that there is a globally asymptotic equilibrium in the positive quadrant and there is no bifurcation.
Let us denote by X the vector field associated to the system (9) and let be the vector field defined on . We denote by the dual form associated with . We can clearly write:
(10) |
(11) |
(12) |
(13) |
(14) |
![](/biologies/article/CRBIOL_2004__327_11_1058_0/jpg/src/tex/figures/gr003.jpg)
Scheme of the return map (Poincaré map) defined for an ɛ-perturbation of a centre. A half straight line is chosen, starting from the equilibrium point. This line is parameterised by a real number h. For each point h on this line, there is a trajectory defined by the system (9), leaving the line and coming back after turning around the equilibrium. The new contact with the line defines the map P(h). If P(h)<h, then the considered trajectory is going to the equilibrium, while if P(h)>h, the trajectory is going away from the equilibrium.
Since for small enough ɛ values and for all h, it results that the equilibrium point is globally asymptotically stable.
5 Conclusion–discussion
A Lotka–Volterra model on two patches has been analysed. One of the patch is a refuge for the prey. This is a way to consider spatial variability. Morevover, the prey growth rates can be different on each patch. In the refuge, the prey has a constant growth rate. It means that if there were no displacement, the prey density would become larger and larger, going to infinity. The first consequence of the displacement is thus a control of the prey density by the predator. Moreover, on the predation patch, the dynamics should be periodic if the displacement were eliminated, since on this patch we considered a Lotka–Volterra model. The equilibrium point should be neutrally stable in that case. We shown that the complete model leads to a globally asymptotic equilibrium. This means that the spatial variability leads to a stabilisation of the dynamics. Even if it is a common result, it is the first time that it is obtained with a mathematical model having appropriate mathematical properties: neutral dynamics on the homogeneous patches and density-independent migration rates. This is the reason why we think that it is a robust result.
More precisely, we can see on Eq. (14) that the larger a, b and are, the more stable the equilibrium is. In other words, when these parameters increase, the trajectories reach the equilibrium faster. a is the predation rate, then a stronger predation leads to a stronger stability. b is the predator growth rate per prey eaten, then a more efficient predator leads to a stronger stability. Finally, is the proportion of prey available for the predators and takes the greatest value when , that is when the prey is uniformally distributed on the patches: a half in the refuge, the other half in the predation patch. This spatial distribution leads to a stronger stability. We can thus conclude that the spatial heterogeneity plays a crucial role in the stabilization, but nor because it could permit the prey to avoid predation, neither because it could enhance the predation. This stabilization is mainly the result of the refuge.
Indeed, if we considered the same kind of model with the predator on both patches (with different rates on each patch), that is without any refuge for the prey, the result would be less obvious. This is the object of a future work.