Revisiting the stability of spatially heterogeneous predator-prey systems under eutrophication

We employ partial integro-differential equations to model trophic interaction in a spatially extended heterogeneous environment. Compared to classical reaction-diffusion models, this framework allows us to more realistically describe the situation where movement of individuals occurs on a faster time scale than the demographic (population) time scale, and we cannot determine population growth based on local density. However, most of the results reported so far for such systems have only been verified numerically and for a particular choice of model functions, which obviously casts doubts about these findings. In this paper, we analyse a class of integro-differential predator-prey models with a highly mobile predator in a heterogeneous environment, and we reveal the main factors stabilizing such systems. In particular, we explore an ecologically relevant case of interactions in a highly eutrophic environment, where the prey carrying capacity can be formally set to 'infinity'. We investigate two main scenarios: (i) the spatial gradient of the growth rate is due to abiotic factors only, and (ii) the local growth rate depends on the global density distribution across the environment (e.g. due to non-local self-shading). For an arbitrary spatial gradient of the prey growth rate, we analytically investigate the possibility of the predator-prey equilibrium in such systems and we explore the conditions of stability of this equilibrium. In particular, we demonstrate that for a Holling type I (linear) functional response, the predator can stabilize the system at low prey density even for an 'unlimited' carrying capacity. We conclude that the interplay between spatial heterogeneity in the prey growth and fast displacement of the predator across the habitat works as an efficient stabilizing mechanism.


Introduction
Understanding top-down control in food webs has been always a focus of theoretical ecology, and revealing potential mechanisms which stabilize predator-prey/resourceconsumer interactions in ecosystems with a high degree of eutrophication (i.e. ecosystems that are high in nutrients) is a particularly challenging problem [42,37,2,22,1,8,31].
Classical theory predicts that a generic predator-prey system will be highly unstable under eutrophic conditions and is expected to exhibit large-amplitude oscillations of species densities which will eventually result in population collapse and extinction [42,19]. However, these theoretical conclusions are often at odds with reality, since there exist a large number of examples of ecosystems in which the species densities remain low despite a high nutrient supply [10,13,7,23,43]. Various mechanisms have been proposed to explain such stability of predator-prey interactions in eutrophic environments. In particular, it has been suggested that this could be due to interplay between structuring within the prey population according to vulnerability or nutrition properties and active food selectivity of the predator (prey switching) [1,17,33,32,31,16]. However, stabilization mechanisms of this kind cease working for unstructured populations, so a natural question is whether efficient top-down control is possible for a population of identical prey individuals.
Spatial heterogeneity of the environment and a non-homogeneous distribution of individuals across the habitat has been suggested to be an important factor promoting top-down control and preventing species extinction [22,36,40,8,39]. In recent model studies of planktonic systems with a high nutrient supply, it was shown that interplay between the gradient of light in the water column and fast vertical displacement of herbivorous zooplankton can stabilize the system even for an unlimited nutrient stock [30,27]. The works cited implement a modelling approach which differs from the classical reaction-diffusion framework in several important ways. In particular, the new approach allows for different time scales for the population growth rate of prey and its predator and also takes into account the fact that the rate at which predator disperses across the habitat can be substantially faster than its characteristic demographic scale. For example, zooplankton herbivores can quickly move vertically along the entire habitat (in the euphotic part of water column) a large number of times between reproduction events, thus we cannot assign an organism to a particular spatial location [11], as in the case of reaction-diffusion models. Thus we should describe the predator population as a whole in terms of the total population size, while still taking into account non-homogeneity of the distribution of predator individuals across the habitat when considering grazing of the prey at each particular location. Ecologically similar scenarios can be found in a large number of other ecosystems when modelling, for instance, interactions between terrestrial herbivorous mammals and grass, zooplankton and fish, insects and birds, insects and their parasitoids, acarine predator-prey systems and even between fish population and fishing boats [29,20,28,24,5,35,3].
Modelling predator-prey interaction in a spatially heterogeneous habitat, where an exact location of individuals on a demographic timescale cannot be properly assigned, requires the implementation of partial integro-differential equations. However, the complexity of this framework makes it rather difficult to treat the model analytically: all previous findings have been obtained by direct numerical simulations of the equations for particular parameterisations of the model ingredients [30,27]. This, obviously, cannot be considered as a rigorous proof of stabilization of the system. Moreover, the results can strongly depend on the choice of the per capita population growth rate, which in the previous studies was considered to be an exponential function modelling light attenuation with depth [30,27]. Finally, we should say that the integro-differential equations in [30,27] differ significantly from those used in the modelling of physiologically structured populations, see e.g. [12,14]. For instance, the model we develop here does not contain a transport term as standard physiologically structured models do.
In this paper, we revisit the previous results on stabilization for a class of integrodifferential predator-prey models, and explore the main factors assuring efficient topdown control. We analytically obtain criteria for the existence of a positive stationary state assuring the coexistence of both species for an arbitrary gradient of resource distribution across the environment. Then we investigate the stability of the coexistence equilibrium and derive a general characteristic equation, which governs the stability of the system. We analytically show that in the case that there is no saturation in the functional response (i.e. Holling type I response), the predator is able to stabilize the system at low prey density even for an 'unlimited' carrying capacity. We study how the size of the habitat and the spatial gradient of the prey growth affects stabilization in the system. We explore two ecologically relevant scenarios: (i) the case when the spatial gradient of the growth rate of the prey is due to abiotic factors only, and (ii) the case where the growth rate of the prey is affected by the surrounding density of the population itself (e.g. effects of self-shading). We find that, interestingly, in case (i) the trophic control of the prey is not possible when the population size of the predator varies slowly, whereas in case (ii) the predator can efficiently suppress the population growth even without changing its own population size.

Modelling framework and biological rationale
The general predator-prey model with a highly mobile predator is described by the following system of partial integro-differential equations [30].
In the model above p(h, t) and z(h, t) denote the densities of prey and predator at time t and at spatial location h, respectively; H is the overall size of the habitat; hence P (t) and Z(t) are the average densities of the species at time t; while p 0 and Z 0 denote the initial conditions. The first term in the right hand side of (2.1) describes the dispersal of prey across the habitat due to diffusion; R is the function (or functional, in general) describing the local per capita growth rate of the prey; f is the local functional response of the predator. The parameters k and m denote the food efficiency of the predator and its mortality rate, respectively. In the case where we model plankton interactions, h is the vertical coordinate from the water surface and H is the overall depth of the euphotic zone.
Here we impose zero flux boundary conditions ∂p ∂h = 0 at h = 0 and h = H. The integro-differential framework (2.1)-(2.3) is quite different from the conventional reaction-diffusion modelling setting. Here we consider that the predator movement occurs on a much faster timescale than both the prey movement and the population dynamics; and assume that the predator can quickly move to any part of the habitat. Thus, we cannot assign a particular spatial location to an individual within the predator population, since organisms constantly change their locations. The framework still allows us to consider the instantaneous spatial distribution of the predator population z(h, t). Note that equations (2.1)-(2.3) do not uniquely specify z(h, t) based on the knowledge of Z(t), so we need to impose some further assumptions. A number of empirical observations suggest that the relative distribution of predator between patches often follows the relative abundance of food, which is described by The above proportion-based parametrization may not be realistic since it does not take into account possible interference between predators [34]. However, in a large number of study cases, see e.g. [29,20,24,26,35], we can still reasonably well approximate the instantaneous spatial distribution of predator across the food patches by (2.4). For instance, the distribution of herbivorous zooplankton in the water column often follows that of phytoplankton, see e.g. [25,30]. Fig.1a shows a typical example of the relative distribution of the biomass of phytoplankton (prey) and herbivorous zooplankton (predator) in different layers of the vertical water column in the euphotic zone of the sea. Here we split the water column into a number of layers and plot the values of the spatially averaged biomass of p and z in each layer divided by the overall average biomass across the entire column. The data were recently collected at 12 stations in the Black Sea for various seasons (2007-2013); further details on obtaining the samples are provided in electronic Appendix 1. The figure presents both night and day observations. For the night observation (open circles) one can see that the data can be roughly described by the straight line with slope approximately equal to 1, thus the basic assumption (2.4) holds (standard linear regression gives R 2 = 0.85 and the slope of the fitting line as 0.96).
The major difference between the night and the day patterns is that the vertical distribution of zooplankton substantially varies due to the phenomenon of regular daily vertical migration: during the day time zooplankton is mainly located in deep layers to avoid predation by higher visual predators (e.g. fish) and ascend to upper layers at night [38]. Note that daytime feeding activity of herbivorous zooplankton is usually low and most of the food is consumed during the night time [6]. Fig.1b provides examples of vertical distributions of species and the daily zooplankton migration pattern. One can see that the vertical distribution of zooplankton is highly variable throughout the 24 hour period. Moreover, zooplankton usually exhibit fast short term unsynchronized migrations [11], by permanently descending and ascending over smaller distances (10-20m). Thus, we cannot assign a particular location to an individual zooplankton and we can use the integro-differential framework with an instantaneous distribution of the predator. We should note that model (2.1)-(2.3) describes the organisms averaged over the day, and therefore excluding daily migration patterns. Finally, Fig.1b also provides the profile of chlorophyll a through the column, which is a proxy of the distribution of phytoplankton. Its complicated structure is due to high heterogeneity in the distribution of the vital resources: light and nutrients (phosphorous and nitrogen) through the column, which potentially translates into a complicated relation r(h) (cf. [41]).
We consider that the local functional response is given by the Holling type II function [21,18] f (p) = α p 1 + pβ , (2.5) where the parameters α and β denote the attack rate and the inverse saturation food density, respectively. In the limiting case when β = 0 the functional response becomes linear or Holling type I response. Here we intentionally implement a concave downwards (i.e. non-sigmoidal) functional response, which is a destabilizing factor in predator-prey models in general, see e.g. [42,19]. For the sake of simplicity we neglect the effect of diffusion of the prey by setting D = 0. In the Discussion section, we shall briefly address the effect that prey diffusion has on the results.
Finally, we need to specify the growth function R. Since we model eutrophic conditions, we consider that this function does not depend on the variation of nutrient concentration. We shall explore two main scenarios: (i) the per capita growth rate is a function of the abiotic environment only, i.e. it is a local function R = r(h) and (ii) the per capita growth rate depends on the population density distribution in some parts of the environment R = R(h, p), i.e. it is a non-local term. In the latter case, we have a functional, i.e. a function which depends on another function characterizing the overall spatial distribution of p. Note that in both cases, we neglect natural mortality of the prey compared to the predation rate.
In case (i) the equations read (Model 1) For the moment we do not specify the function r. We only assume that it is a smooth and non-negative function.
In case (ii), we shall investigate the particular scenario of self-shading in the population. In this complex scenario, the supply of a resource such as light decreases in space via a function of the local population density. In planktonic communities this may be shading of light for phytoplankton residing in deeper layers by the organisms dwelling in upper layers. For a non-planktonic system, this may be for instance a gradual depletion of a nutrient as it flows through part of the environment occupied by the prey and is gradually consumed by the organisms. In this case, the local rate of decrease of resource density can be assumed to be a function of the density of consumers. It is natural to assume that the available resource quantity is a decreasing function of the total number of consumers h 0 p(x, t) dx that are above/upstream an individual at depth/point h.
In particular, we assume that the quantity of the available resource is an exponentially decreasing function of the number of consumers that contribute to the shading.
The equations read (Model 2) where ν is the coefficient of self-shading. In both of the models above the average densities P (t) and Z(t) are given by equation (2.3).
In this paper, we are mostly interested in the existence and stability of the positive stationary state of (2.6)-(2.7) and (2.8)-(2.9). It is worth mentioning that a spatially homogeneous predator-prey model (2.1)-(2.3) with perfect mixing (i.e. for r(h) = const and z(h) = const) will be globally unstable [4], thus the main question is whether it is possible to stabilize this system by introducing spatial heterogeneity (i.e. structuring the population) and allowing a non-homogeneous distribution of predators. It is also worth noting that the term describing prey growth in equation (2.8) has infinite dimensional nonlinearity, and that structured population models with such nonlinearity are notoriously difficult to analyse, see e.g. [9], [15].

Model 1:
The local growth rate of the prey is a function of abiotic factors only.

Existence of a non-trivial (coexistence) steady state.
We look for a strictly positive steady state (p * , Z * ) of model (2.6)-(2.7). We assume that r is a strictly positive continuous function with We characterise the existence of the positive steady state in the following proposition. Proof. The steady state equations read For positive values (P * , Z * ) which satisfy we obtain from equation (3.11) We then substitute the steady state p * given by formula (3.15) into equations (3.12)-(3.13). After simplification we obtain: From equation (3.17) we obtain: from which we obtain We substitute this expression into equation (3.16) to obtain Also note that (3.14) together with (3.19) require the existence of a (unique) solution of (3.20) in the interval Hence there is a solution of equation (3.20) in the required interval (3.21). It is also shown that on the interval (3.21) the function F is strictly monotonically increasing, hence the solution is unique, and the proof is complete. Note that, we may relax the assumption of strict positivity of r, and allow it to vanish at some values of h. Then equation (3.15) shows that the steady state p * will also vanish at those values of h, but Proposition 3.1 still holds.

Stability of the coexistence steady state.
The stability of the coexistence equilibrium is determined by the roots of a characteristic equation, which we derive in the current subsection. The linearisation around the coexistence steady state (p * , Z * ) of model (2.6)-(2.7) is computed as We look for the solution in the form of ∆p(h, t) = w(h) exp(λt) and ∆Z(t) = U exp(λt). This leads to the eigenvalue problem (W is the spatial average of w) Using (3.15) we recast equation (3.25) as follows We re-arrange equation (3.27) as Note that k 1 (λ, ·) = 0 for λ ∈ C, Re(λ) ≥ 0. Next we integrate equation (3.30) from 0 to H and divide by H to obtain where We also multiply equation (3.30) by r(h) H and integrate from 0 to H to obtain where dh.
Next we substitute into equation (3.28) the expression for W given by (3.32). This, together with equation (3.31), yields a homogeneous system of two equations for the variables (W, U ) as follows.
For any value of λ ∈ C with Re(λ) ≥ 0 a non-trivial solution of system (3.33) yields a non-trivial solution w via equation (3.30). On the other hand, any eigenvalue λ with Re(λ) ≥ 0 necessarily satisfies equations (3.33) for same non-trivial (W, U ). We summarize our result in the following theorem: Theorem 3.2. The non-trivial steady state (p * , Z * ) of model (2.6)-(2.7) is locally asymptotically stable if every solution λ of the characteristic equation below has negative real part.
where k 2 and k 3 are defined earlier. On the other hand, the steady state is unstable if there is at least one solution of (3.34) with positive real part.
In general, the derived characteristic equation (3.34) is complicated and determining its roots might require extensive numerical simulations which might take as much time as directly solving the system of differential equations. However, this equation can be rather useful for construction of a Hopf bifurcation curve: one can set λ = ib (b is a positive real number) and solve the resulting equation. In the Discussion section we shall provide an example of construction of such a curve for a particular parametrization of the function r.
Using the characteristic equation above we are able to address the stability of the system in the case when the demographic timescale of the predator is much slower than that of the prey. In this case, we can assume the population size of the predator to be constant. The result is given by the following remark.
It is shown that k 2 (0) > 1, and clearly k 2 is a monotone decreasing function of λ for λ > 0. Also we have that lim λ→+∞ k 2 (λ) = 0, hence there exists a unique positive eigenvalue, and the steady state is unstable.
Remark 3.4 Another important conclusion we may draw from (3.34) is that for a homogeneous environment (i.e. r(h) ≡ const) the system is always unstable provided β > 0, since all eigenvalues will have positive real parts (we do not show here the corresponding characteristic equation for the sake of simplicity).

3.1.3.
Stability in the special case β = 0. Next we investigate an important special case of (3.34) where we have a Holling type I functional response, i.e. there is no saturation rate in predation, β = 0. The characteristic equation reduces to which can be simplified by recalling the expressions for the stationary densities of species (3.16), (3.17) and (3.19). We obtain Let us first consider the case of a small environment, i.e. small H. We use the Taylor expansion of (3.37). We assume that the derivative of r(h) at h = 0 does not vanish: r(h) ≈ r 0 + Rh, where R = 0. The characteristic equation becomes (3.38) This expression can be re-written in the polynomial form (λ = −r 0 ) by dropping the We use the Routh -Hurwitz stability criterion for the coefficients of the above polynomial: a n > 0, a 2 a 3 > a 4 a 1 , a 1 a 2 a 3 > a 4 a 2 1 + a 0 a 2 3 . One can see that for a sufficiently small H: a n > 0 (n = 1, 2, 3, 4), moreover Thus, for small H the polynomial characteristic equation (3.39) has eigenvalues with only negative real parts and, consequently, the initial characteristic equation (3.37) also exhibits stability for small H. Note that stability is sustained in the case when we consider a parabolic approximation r(h) ≈ r 0 + Rh + R 1 h 2 , where R = 0, R 1 = 0 (the need for considering such an approximation comes from the fact that in the second order expansion (3.38) the quadratic term R 1 will be present. We should stress that here we are allowed to drop O(H 3 ) and consider (3.39) since all the roots of (3.39) are perturbations of the solutions with H = 0 such that the maximal order of perturbation with respect to H does not exceed 2, i.e. it is less than 3 (we do not show the corresponding expressions λ(H) of solutions of (3.39) for the sake of brevity). Now let us suggest that for some critical H > 0 the stability of (3.37) changes and the real part of an eigenvalue becomes positive. This can occur at a Hopf bifurcation point and the condition of stability change is λ = ib, where b is a positive real number. We substitute λ = ib into the complex characteristic equation and obtain two equations for its real and imaginary parts. After some re-arrangement the equation for b becomes (see electronic Appendix 2 for the exact derivation) We substitute a particular parametrization r(h) to verify whether or not (3.42) has a real value solution. In the case (3.42) does not have real solutions, the steady state of (2.6)-(2.7) will be always locally stable for β = 0 because for small H the eigenvalues always have non-zero negative real parts and for any H > 0 they will not vanish (otherwise we would have a solution with λ = ib, in which case b would solve (3.42)). However, it is still an open question whether or not the equilibrium will remain stable for β = 0 for any arbitrary positive function r(h). In our investigation, we checked several concrete families of ecologically relevant and mathematically tractable functions including linear, parabolic, cubic, exponential and sinusoidal functions given, respectively, by r(h) = a 0 + a 1 h, r(h) = a 0 + a 2 h 2 , r(h) = a 0 + a 3 h 3 , r(h) = a 0 exp(a 4 h) and r(h) = a 0 (1 + sin(ωh)), where a i , ω, are parameters. For those functions we can obtain analytical expressions for the corresponding integral and numerically check the possibility to find a solution of (3.42), by varying the parameters. For all these functions we find that (3.42) has no real solution, thus the equilibrium is always stable (we do not show the corresponding cumbersome expressions and graphs here for the sake of brevity).
3.2. Model 2: local growth rate of prey is affected by self-shading.
We look for the coexistence equilibrium (p * , Z * ) of model (2.8)-(2.9). The steady state equations read From the first equation of the system we obtain . (3.46) We differentiate the above expression to obtain Integration of (3.47) gives where C is a constant of integration. The constant C can be found by setting h = 0 and using the expression for p * (0) from (3.43) p * (0) = P * r 0 αZ * − r 0 P * β , (3.49) After substituting C and some simplification, we obtain the expression for the stationary state p * (h): The above expression is implicit since it is impossible to solve it analytically for p * (h) except in the particular case β = 0; however, from (3.47) one can conclude that p * (h) is a decreasing function of h. The maximal density through the column is attained at h = 0 and given by (3.49); the minimal density is achieved at h = H. From (3.43) we obtain the expression for p * (H) We use (3.51) and (3.50) to get the following relation between P * and Z * 0 = − exp(P * Hν)αZ * + βP 2 * Hν + Hνr 0 P * + αZ * . (3.52) Next we consider the stationary equation (3.44) for the predator and replace the integration variable h with p(h) using (3.47).
We are allowed to make the above change of variables since p(h) is a monotonically decreasing function. We substitute the values of p * (0) and p * (H) from (3.49) and (3.51), respectively, to obtain after simplification another identity linking P * and Z * Z * = (1 − exp(−P * Hν))r 0 k mHν . It is easy to see that (3.56) will be satisfied for some P * > 0 if and only if αk − βm > 0. Further, by computing the second derivative of f (P * ) one can show that for αk − βm > 0 it is always negative, thus there can be only one positive solution of (3.55). Hence, combining the above properties of f (P * ) and f 1 (P * ), we find that the existence of the nontrivial steady state of the system (2.8)-(2.9) is given by the following It is easy to see that for β 1 the conditions of the proposition are satisfied, thus there always exists a non-trivial steady state. Let us consider a gradual increase in β while still keeping αk −βm > 0. The root of f (P * ) = 0 will increase, whereas the root of f 1 (P * ) = 0 will decrease. At a certain β * both roots coincide, which is an unavoidable event since the root of f 1 (P * ) = 0 will approach zero when αk = βm. Thus for all β < β * , the system has a non-trivial coexistence equilibrium whereas for β > β * , a positive equilibrium is impossible. The exact value of β * can be found from the condition f (P * ) = f 1 (P * ) = 0 which results in a cumbersome analytical expression (not shown here).

3.2.2.
Stability of the non-trivial (coexistence) steady state of Model 2.
The linearisation around the positive steady state (p * , Z * ) of (2.8)-(2.9) results in a rather cumbersome characteristic equation. In this paper we focus on an important particular case, where there is no saturation in the functional response of predator, i.e. β = 0. The linearized system becomes We look for the solution in the form of ∆p(h, t) = w(h) exp(λt) and ∆Z(t) = U exp(λt). We substitute the expressions for the stationary solutions (3.48) (β = 0) and (3.54) to arrive at the following eigenvalue problem (W is the spatial average of w) We re-arrange the first equation of the above system We differentiate both sides of (3.61) and obtain dw dh After simplifying the above equation we get We suggest that (hν + C)λ + νC = 0 and derive the following equation for the eigenfunction w(h) Integration of this linear equation gives the following eigenfunction (3.65) Next, we substitute this function into system (3.59)-(3.60) where the value of W is obtained from . (3.68) We simplify equation (3.67) finally, we substitute the expressions for the stationary value P * into (3.66) and (3.69) to obtain the characteristic equation where Q 1 = Hν (C + Hν) ln((C + Hν)/C) < 1, (3.71) It is easy to see that characteristic equation (3.70) always has negative real parts, thus, the linearized system (3.57)-(3.58) is always stable. We should also say that unlike the model without self-saturation the number of eigenvalues is finite and there is a single eigenfunction. This can be explained by the fact that separation of variables provides us with only the solution for large values of time and does not describe transient regimes at small and intermediate times.
The above result can be summarized as the following Theorem 3.6. The steady state (p * , Z * ) of model (2.8)-(2.9) is always locally asymptotically stable for β = 0 (Holling type I functional response).
Remark 3.7 One can prove that in a similar way the stability of equation (3.57) governing the dynamics of prey for a constant (equilibrium) density of predator, i.e. for Z(t) = const or U (t) = 0. The eigenvalue of the system is given by This signifies that a spatial regrouping of the predator population (without variation in the total population size) alone can control the prey population in the case of selfshading.

Discussion and summary of results
In this paper, we revisit predator-prey interaction in ecosystems with a high nutrient load and explore the role that spatial heterogeneity of the environment and high predator mobility have in stability of such ecosystems. Our analytical investigations of the model predict that interplay between these two factors results in stabilization of an otherwise globally unstable predator-prey system. An important assumption for the realization of the given stabilization scenario is that the predator preferentially feeds at locations with high prey density (see relationship (2.4)). Since spatial gradients in a species' growth rate across their habitats is a typical ecological situation, and many predators/consumers are highly mobile and feed mainly in patches of high food density [29,20,28,24,5,35] (see also figure 1 of the current paper), we claim here that the reported mechanism should be rather typical for marine, freshwater and terrestrial ecosystems under eutrophication. Our findings also provide a new and robust mechanism to resolve the famous 'paradox of enrichment', which is a long standing enigma in theoretical ecology [42,19,1,17,39,33,43].
Our models are formulated as integro-differential equations, which allows for the incorporation of different time scales: fast spatial movement of predator and slow population dynamics of both species. This makes the current framework more effective than classical reaction-diffusion models for modelling trophic interaction in small-sized habitats, since individuals cannot be assigned to a particular spatial location in this case while reaction-diffusion equations contain only local terms.
Here we explore a generic predator-prey model with an infinitely large carrying capacity of prey (i.e. no resource limitation) and consider two main scenarios: (i) the prey growth rate gradient is fixed (Model 1) and the local growth rate of prey is affected by self-shading (Model 2). The main mathematical outcomes of our study are the following: (i) For each model (Model 1 and Model 2) we find the condition for the existence of the positive stationary state through which prey and predator can coexist for an arbitrary gradient of the growth rate r(h) ≥ 0. Interestingly, including self-shading usually restricts the conditions of coexistence (cf. Propositions 3.1 and 3.5).
(ii) For Model 1 we have derived the general characteristic equation (3.34) governing the stability of the coexistence equilibrium. Using this equation, one can construct the Hopf bifurcation curve that separates stable and unstable dynamics without numerically integrating the underlying partial integro-differential equations (see an example below).
(iii) Analysis of the characteristic equation (3.34) shows that for a homogeneous environment, in which r(h) ≡ const, the system is always unstable provided there is saturation in the functional response (Remark 3.4). Thus, spatial heterogeneity of the environment is a necessary condition for stability of the coexistence equilibrium.
(iv) For both Model 1 and Model 2, in the case of a linear (Holling type I) functional response of the predator, we analytically proved that the coexistence equilibrium is always stable. Since the system is structurally stable, the stability of the equilibrium is also guaranteed for sufficiently small β > 0, i.e. for a Holling type II functional response with low saturation.
(v) Comparison of Model 1 and Model 2 shows that stabilization in the system with dynamical self-shading (Model 2) is potentially a stronger mechanism than the scenario, where the spatial gradient in the prey growth rate is fixed (Model 1). For instance, without saturation in the functional response (β = 0) the stabilization does not require changing the predator biomass: variation in the spatial distribution of a fixed number of predators can suppress the prey outbreak and re-establish the equilibrium (cf. Remarks 3.3 and 3.7).
In the case of a fixed prey growth rate gradient, the general characteristic equation (3.34) allows us to easily compute the Hopf bifurcation curve to direct numerical simulation of the underlaying integro-differential equations (2.6)-(2.7). As an illustrative example, here we have constructed the Hopf bifurcation curve for the system with a linear variation of the growth rate r(h) = a 0 + a 1 h (see electronic Appendix 3 for detail). Fig.2 shows a set of bifurcation curves in the diagram a 1 and β constructed for a gradually increasing size H of the habitat. In each case, the stability region is located below the curve. One can see that increased the saturation in food consumption, signified by a large β, usually acts as a destabilizing factor, but the influence of the gradient of the growth rate a 1 is more complicated: small and large a 1 result in destabilization of the system whereas stabilization occurs for some intermediate a 1 . An increase in the size of the habitat can also either destabilize or stabilize the system. This result is somewhat counterintuitive since in earlier works it was shown that for an exponential parametrization of r(h), an increase in the spatial gradient of the growth function or the size of the habitat will always facilitate stabilization of the system [30,27].
Our analytical investigation confirms the numerical results reported earlier in [30,27] about the possibility of the stabilization of systems with 'unlimited carrying capacity', and extend the previous findings to more general parametrizations of the spatial gradient of r(h). In particular, we proved that for a small habitat size H, any arbitrary nonhomogeneous dependence r(h) will stabilize the predator-prey interactions (see Section 3.1.3). On the other hand, our work shows that using some other parametrizations, for example a linear function, can provide some qualitatively different prediction compared to the exponential dependence considered in [30,27]. This underlines the danger of sticking to particular parametrizations when modelling predator-prey systems.
Throughout the paper we have neglected the diffusion of the prey population by setting D = 0 in equation (2.1). Including diffusion D > 0 will make the model more realistic, but the equations become analytically untractable. Although one can still perform linearization techniques, it becomes impossible to obtain the resultant characteristic equation due to the presence of the second order derivative with respect to x, thus the only way of exploring the system is by direct numerical simulation. For small values of D the main conclusions of our paper obtained for D = 0 will hold due to the structural stability of the model equations. Interestingly, for some ecosystems (e.g. planktonic systems, [30,27]) the diffusion coefficient D > 0 is relatively small (compared to the interaction terms) so the results obtained should be close to those with D = 0 . However, for a very large diffusion coefficient corresponding to a highly mobile prey, the system becomes homogeneous in space and is reduced to the classical Rosenzweig-MacArthur model, which is always globally unstable for any β > 0. Thus, the system will be destabilized for a certain critical D, which can be found numerically.
Finally, we should mentioned the open questions that our study raises. For instance, it would be an important extension to determine the stability or otherwise of the coexistence equilibrium of Model 1 in the case r(h) is an arbitrary function for a Holling type I functional response (see equation (3.42)). Based on preliminary results obtained for a small habitat size H, we hypothesize that for β = 0 the Model 1 should be stable for an arbitrary r(h), but it remains to prove this conjecture rigorously. Another challenging issue is exploring the system stability for a more complicated relationship between z(h, t) and p(h, t), by taking into account possible interference between predators [34]. For instance, this can be represented by setting z(h) = p(h) µ P µ Z, µ > 0 [34]. Another interesting extension would be to investigate analytically the model combining self-shading of the population with a fixed gradient of local growth rate and combining Model 1 and Model 2. From the biological point of view it might even be natural to consider situations when the growth function R (or r) can take negative values, i.e. the population of prey actually declines locally. In this more general case, model (2.6)-(2.7) will still be well-posed, and solutions starting in the positive cone would remain non-negative for all times. At the same time equation (3.11) shows that there will be no strictly positive steady states of the model and this might influence the previous stability results. Models with locally negative growth rate of prey could therefore possess quite rich dynamics, and such models should definitely be a matter of further investigation.