Introduction

Several predator–prey species (Murdoch et al. 2002) show oscillatory dynamics and some predator–prey cycles exhibit different amplitudes across their geographic range (Ruggiero et al. 2000). It has been suggested that geographic differences in cycle amplitude may be due to variations in the extent of habitat fragmentation from one region to another (Ruggiero et al. 2000; Ylonen et al. 2003). Habitat fragmentation produces spatial heterogeneity in the domain, and the importance of spatial heterogeneity in ecology has been recognized for the past few decades (Kareiva 1990). There are many different concepts which have linked spatially non-uniform habitats to population dynamics, such as metapopulations (Levin 1974), island–mainland populations (see references in (Kareiva 1990)), the refugium hypothesis (Ruggiero et al. 2000), and source-sink populations (Pulliam 1988). Empirical studies have also found connections between habitat fragmentation and the dynamics of mammalian cycles, such as snowshoe hares (Akcakaya 1992; Murray 2000) and meadow voles (Ylonen et al. 2003).

The effects of habitat fragmentation on cyclic populations have been reviewed by Fahrig 2003 and Ryall and Fahrig 2006. In their reviews, Fahrig and Ryall observed that habitat can be fragmented in a number of different ways (Fahrig 2003) and that these different fragmentation types can vary in their effect on species persistence. Fragmented habitat is characterised by the relative size of high- (good) and low-quality (bad) habitat patches and the distances between them. Theoretically, fragmentation can be increased by decreasing good patch size, increasing bad patch size, or by making both changes simultaneously. By categorizing fragmentation research by fragmentation type, Fahrig and Ryall were able to discern trends in the effects of habitat fragmentation. They found that decreasing good patch size had negative effects on biodiversity (Fahrig 2003) and had a greater negative effect on the abundance of a specialist predator than on their prey (Ryall and Fahrig 2006). These trends in population dynamics in response to habitat fragmentation were consistent with the results of Strohm and Tyson 2009.

Strohm and Tyson 2009 used diffusion–advection–reaction models to investigate the effects of habitat fragmentation on cyclic population dynamics. They found that increased habitat fragmentation results in a decrease in cycle amplitude and changes in average population densities. In particular, habitat fragmentation results in a decrease in the average density of predators while the average density of prey may increase or decrease depending upon the choice of the reaction terms.

These results were chiefly numerical, since the spatially explicit models were not amenable to analytic approaches. There have been techniques developed to reduce spatially explicit models into spatially implicit equivalents (Dieckmann et al. 2004; Tilman and Kareiva 1997). These techniques include the aggregation approach, average dispersal success, and the use of modified mean-field equations (Bascompte 2001; Cobbold 2005). In one particular study, Cobbold et al. 2005 found that they could successfully reproduce the critical patch size dynamics and the bifurcation structure of a spatially explicit model with a spatially implicit model using average dispersal success (Van Kirk and Lewis 1997).

The goal of this study is to find a simple spatially implicit ordinary differential equation (ODE) model that can explain the response of the spatially explicit partial differential equation (PDE) model (Strohm and Tyson 2009) to habitat fragmentation. We develop the spatially implicit models by first formulating a related but simpler PDE problem on a single patch with Dirichlet boundary conditions and then by extracting the principal eigenvalue derived from a first-term Fourier truncation of the simpler PDE model. The principal eigenvalue of a spatially explicit model is a measure of persistence of the population (Cantrell and Cosner 2003). The simplified ODE models are analysed and provide some insight into the behaviour of the PDE models. We find that habitat fragmentation produces similar responses in both models for larger high-quality patch sizes. Furthermore, the critical patch size for both predator and prey and the densities at which the Hopf bifurcations occur in the ODE models provide an upper bound for equivalent behaviour in the PDE models. We discuss the relationship between the spatially implicit ODE model and original three-patch PDE model and show that the differences in behaviour can be anticipated based on the nature of our simplification procedure. In particular, we confirm that the introduction of Dirichlet boundary conditions on a single patch mimics lethal boundaries which we expect to lead to an overestimate of critical patch sizes in the PDE model.

In “Spatially explicit models”, we review the modelling framework for the spatially explicit models. Then in “Simplification to ODEs”, we introduce our simplification technique which is based on the method of Fourier series (Keener 2000) and reduce the spatially explicit models to spatially implicit ODE models. The spatially implicit and spatially explicit models are solved, using the methods in “Methods”, and the results of these numerical simulations are compared in “Comparison of ODE and PDE results”. “Analysis of the ODE modelsAnalysis of the ODE models” provides analytical results on the critical patch size and steady states in the ODE models. This section also includes bifurcation analysis and nullcline plots to analyse the behaviour of the coexistence steady states in the spatially implicit models. Finally, our results and future directions for research are summarized in “Discussion”.

Spatially explicit models

The spatially explicit PDE models are taken from Strohm and Tyson 2009 and take the form

$$ \frac{\partial n}{\partial t} = D_n \frac{\partial^2 n}{\partial x^2} + f(n,p,x), \label{eq:primary-a} $$
(1a)
$$ \frac{\partial p}{\partial t} = D_p \frac{\partial^2 p}{\partial x^2} + g(n,p), \label{eq:primary-b} $$
(1b)

where,

D n , D p :

The diffusivity coefficients for prey and predator, respectively,

n, p :

The population densities of prey and predator, respectively,

f(n,p,x), g(n,p):

The reaction terms for prey and predator, respectively.

Note that the reaction terms f(n,p,x) in Eq. 1a depend explicitly on position, x. This is because we allow the prey growth rate to vary spatially with habitat quality. Specifically, the prey growth rate is zero in bad patches and a positive constant in the good patches. To ensure that r(x) is continuous, we smoothly decrease the growth rate to zero at the edges of the good patch.

We take the functions f(n,p,x) and g(n,p) from the historical literature of predator–prey models (Strohm and Tyson 2009). These are the Lotka–Volterra (LV) (Murray 1993; Okubo and Levin 2001), Rosenzweig–MacArthur (RM) (Turchin 2001; Rosenzweig and MacArthur 1963), May (May 1974), and variable territory (VT) (Turchin and Batzli 2001) models. In this study, we depart somewhat from the models studied by Strohm and Tyson 2009 in that we exclude advection-type terms (for simplicity) and use reaction terms with an alternate form for logistic growth (Cosner 1996; Enright 1976) in the RM, May, and VT models. The logistic term used in Strohm and Tyson 2009 is written

$$ r(x)n\left( 1-n/k \right), \label{eq:past} $$
(2)

we instead use

$$ n\left(r(x)-n/k \right). \label{eq:alternate} $$
(3)

This change makes little difference to the PDE model results but is important in the spatially implicit models we develop here (“Spatially implicit models”). We discuss the reasons for using the alternate form (Eq. 3) in “Spatially implicit models”. Mathematically, we have

  1. LV
    $$ f(n,p,x) = r(x)n-cnp,\label{eq:lv-a} $$
    (4a)
    $$ g(n,p) = \chi cnp-\delta p,\label{eq:lv-b} $$
    (4b)
  2. May
    $$ f(n,p,x) = n\left( r(x)-\frac{n}{k}\right) -\frac{cnp}{d+n},\label{eq:may-a} $$
    (5a)
    $$ g(n,p) = p\left( s-\frac{qp}{n}\right),\label{eq:may-b} $$
    (5b)
  3. RM
    $$ f(n,p,x) = n\left( r(x)-\frac{n}{k}\right) -\frac{cnp}{d+n},\label{eq:rm-a} $$
    (6a)
    $$ g(n,p) = \frac{\chi cnp}{d+n}-\delta p,\label{eq:rm-b} $$
    (6b)
  4. VT
    $$ f(n,p,x) = n\left( r(x)-\frac{n}{k}\right) -\frac{cnp}{d+n},\label{eq:vt-a} $$
    (7a)
    $$ g(n,p) = \frac{\chi cnp}{d+n}-\delta p -\frac{sqp^{2}}{n},\label{eq:vt-b} $$
    (7b)

Parameter values

For comparison of the model solutions, we need to choose a set of parameter values for which the models, Eqs. 1a and 1b with Eqs. 4a–7b, yield cycles with comparable periods and prey average values. We therefore select the parameter values used by Strohm and Tyson 2009. These were motivated by the Canada lynx and snowshoe hare predator–prey system in the boreal forest. Parameter values are listed in Table 1.

Table 1 Parameters and default values for the LV, RM, May, and VT reaction terms
Full size table

Simplification to ODEs

Strohm and Tyson (Strohm and Tyson 2009) numerically investigated and compared the solution behaviour of the four PDE models, Eqs. 1a and 1b with Eqs. 4a–7b. Analytic solutions were not available, and so the mechanism behind the observed trends in population cycle amplitude and average densities was not evident. Below we describe a method for simplifying the PDE system to an analytically tractable ODE one, and compare the solutions from the two classes of models. The parameters used in the ODE models are the same as those used in the PDE models. To perform the simplification, we assume that the multi-patch spatially explicit PDE model can be approximated by a single-patch spatially explicit model with homogeneous Dirichlet boundary conditions. In Appendix B, we also investigate the single-patch spatially explicit model with partially absorbing boundary conditions, which more closely approximate the original three-patch problem.

In our equations (Eqs. 1a and 1b), we consider movement governed by diffusion only, and so the spatial operators are self-adjoint. Thus, the principal eigenvalue for the solutions can be calculated using Fourier series (Kielhofer 2004). Furthermore, we note that the bifurcation structure of the nonlinear system (Eqs. 1a and 1b) is exactly determined by replacing the spatial operators in Eqs. 1 and 2 with their first eigenvalues (Keener 2000; Kielhofer 2004), obtained by linearizing about appropriate steady-state solutions. This is the general framework of our approach. We will find that linearization about the nontrivial steady states does not lead to analytic insights but, surprisingly, linearization about the extinction steady state is highly informative even at states well removed from the extinction state. Linearization about the nontrivial steady states of the system and principal eigenvalue analysis is a common tool in investigating persistence of predator and prey (Cantrell and Cosner 2003). By investigating the principal eigenvalue of the spatially explicit models, we are able to determine the appropriate manner to include the diffusion into the growth/death rate of the corresponding spatially implicit model.

Model simplification

To obtain a spatially implicit model, we first consider the most basic spatially explicit models, Eqs. 1a and 1b with Eqs. 4a and 4b, on a single good patch with r(x) =constant. If we consider the prey in the absence of predators, we have

$$ \frac{\partial n}{\partial t} = rn+D_n\frac{\partial^2 n}{\partial x^2}. \label{eq:simple} $$
(8)

We begin our analysis by considering Eq. 8 on an interval 0 ≤ x ≤ L with homogeneous Dirichlet boundary conditions,

$$ n(0,t)=0, \quad n(L,t)=0. \label{eq:bc} $$
(9)

The same equation (Eq. 8) is obtained when using RM, May, or VT reaction terms by linearizing around the trivial steady state. The solution to Eqs. 8 and 9 is found by the method of Fourier series (Haberman 2004). This provides the solution

$$ n(x,t)=\sum\limits_{m=1}^{\infty} A_m e^{[r-D_n(m \pi/L)^2]t}\sin\left(\frac{m\pi x}{L}\right), \label{eq:eigfn} $$
(10)

with eigenvalues \(\lambda_m =r-D_n(m\pi/L)^2\). The m = 1 term dominates all other terms in this infinite series and so we drop all higher order terms. The λ 1 parameter is known as the principal eigenvalue, and can be used to determine the persistence of the population (Cantrell and Cosner 2003): if λ 1 > 0 the population grows exponentially (persistence) while if λ 1 < 0 the population decays exponentially (extinction).

We are interested in the cases leading to population persistence, and therefore, we consider only the solution term in Eq. 10 corresponding to the first eigenvalue, that is

$$ A_1 e^{[r-D_n(\pi/L)^2]t}. $$
(11)

This quantity is the exact solution to the ODE

$$ \frac{{\rm d}n}{{\rm d}t} = \left(r-D_n\left(\frac{\pi}{L}\right)^2\right)n. \label{eq:newode} $$
(12)

We thus arrive at the result that the leading-order population dynamics are the solution to Eq. 12 with growth rate \(r-D_n(\pi/L)^2\), which varies with diffusion rate, D n , and the domain size, L. In Eq. 12, the growth rate of the population will go to zero as L decreases, giving rise to a critical patch size, L c , below which the population cannot persist. The growth rate of Eq. 12 is the same as the principal eigenvalue \(\lambda_1 = r-D_n(\pi/L)^2\) of the spatially explicit model with Dirichlet boundary conditions, and so λ 1 can be used to determine population persistence in the spatially explicit model (8) with (9).

In the absence of predation, we have steady-state prey populations given by the equilibrium solutions to Eq. 8. We refer to these as \(\bar{n}(x)\). In order to determine if the predator population can persist, the usual procedure is to linearize about the zero predator and prey equilibrium solution \(\bar{n}(x)\) and determine the subsequent growth rate of the predator population. We carry out this approach in Appendix A and show that the analytic solution cannot be obtained. So we take an alternate approach, where we linearize the model about (n,p) = (n  ∗ ,0) where n  ∗  is a spatially unvarying prey population density. This gives us

$$ \frac{\partial p}{\partial t} = \chi c{n^\ast}p-\delta p +D_p\frac{\partial^2 n}{\partial x^2}. \label{eq:simplep} $$
(13)

As before, we consider Eq. 13 on an interval 0 ≤ x ≤ L with homogeneous Dirichlet boundary conditions,

$$ p(0,t)=0, \quad p(L,t)=0. \label{eq:bcp} $$
(14)

In this way, we find that the leading-order population dynamics of the predator is governed by the ODE

$$ \frac{{\rm d}p}{{\rm d}t} = \left(\chi c {n^\ast}-\delta -D_p\left(\frac{\pi}{L}\right)^2\right)p. \label{eq:newode2} $$
(15)

In the case of Eqs. 1a and 1b with RM and VT reaction terms, we arrive at a similar expression with χc n  ∗  replaced by \(\frac{\chi c{n^\ast}}{d+{n^\ast}}\).

Spatially implicit models

Based on the result (Eq. 12), we assume that the growth rate of the prey equation in all four models varies with D n and L in the same way. Similarly, we assume that the result (Eq. 15) applies to the death rate in the predator equation for the LV, RM, and VT equations. We define

$$ r_L = r-D_n\left( \frac{\pi}{L} \right)^2. \label{eq:r} $$
(16a)
$$ \delta_L = \delta+ D_p \left( \frac{\pi}{L} \right)^2. \label{eq:delta} $$
(16b)

Since the May equation does not possess a predator death rate, it is the predator growth rate that varies with D p and L. Linearizing about (n  ∗ ,0) as above, we find

$$ s_L = s- D_p \left( \frac{\pi}{L} \right)^2. \label{eq:s} $$
(17)

We thus arrive at the following set of spatially implicit models

  1. LV
    $$ \frac{{\rm d}n}{{\rm d}t} = r_Ln-cnp,\label{eq:lvode-a} $$
    (18a)
    $$ \frac{{\rm d}p}{{\rm d}t} = \chi cnp-\delta_Lp.\label{eq:lvode-b} $$
    (18b)
  2. May
    $$ \frac{{\rm d}n}{{\rm d}t} = \left( r_L-\frac{n}{k}\right)n -\frac{cnp}{d+n}, \label{eq:mayode-a} $$
    (19a)
    $$ \frac{{\rm d}p}{{\rm d}t} = p\left(s_L-\frac{qp}{n}\right).\label{eq:mayode-b} $$
    (19b)
  3. RM
    $$ \frac{{\rm d}n}{{\rm d}t} = \left( r_L-\frac{n}{k}\right)n -\frac{cnp}{d+n},\label{eq:rmode-a} $$
    (20a)
    $$ \frac{{\rm d}p}{{\rm d}t} = \frac{\chi cnp}{d+n}-\delta_Lp.\label{eq:rmode-b} $$
    (20b)
  4. VT
    $$ \frac{{\rm d}n}{{\rm d}t} = \left( r_L-\frac{n}{k}\right)n -\frac{cnp}{d+n},\label{eq:vtode-a} $$
    (21a)
    $$ \frac{{\rm d}p}{{\rm d}t} = \frac{\chi cnp}{d+n}-\delta_Lp -\frac{sqp^{2}}{n}.\label{eq:vtode-b} $$
    (21b)

These spatially implicit models govern the exact bifurcation structure of Eqs. 1a and 1b on a single good patch of length L with homogeneous Dirichlet boundary conditions (Kielhofer 2004). We will compare the solution behaviour of these models with the original three-patch PDE model (“Comparison of ODE and PDE results”).

Recall that the logistic growth terms in the models above are written in an alternate form where the growth rate r L appears inside the density-dependent term. In the spatially explicit models, this growth rate was r(x), whereas here we have r L , which depends explicitly on L but not on x. Thus, the carrying capacity of prey (and predators in the spatially implicit May model) changes with patch size L. Biologically, this means that smaller patch sizes will support a lower density of individuals at a given point inside the domain. In contrast, if the standard form of the logistic equation (Eq. 2) is used, the carrying capacity does not change with patch size, L. This behaviour is problematic since it implies that the carrying capacity for each population is independent of L. For this reason we adopt the alternate form of the logistic equation for this study.

Note that we have derived the functional dependence of the growth and death rates on patch size L for the case of a single patch with Dirichlet boundary conditions. We are using this problem to approximate the original three-patch domain consisting of a good patch surrounded by two bad patches. In the original problem, the bad patches are areas in which there is simply no growth of prey: death rates are not increased in bad patches. We thus expect that the population densities at the boundaries between the good patch and surrounding bad patches to be positive, reflecting some movement of prey and predators between the different patches. In our simplification of Eq. 8 with Eq. 9, the Dirichlet boundaries on the single good patch correspond to lethal boundaries between the good and bad patches, that is, the prey and predator populations immediately drop to zero at the boundaries. Consequently, we expect that the critical patch size predicted by analysis of the spatially implicit models with Eqs. 16a–17 will be an overestimate of the critical patch size of the corresponding spatially explicit models on the original three-patch domain. This suggests that partially absorbing (Robin) boundary conditions on a single patch might give a better approximation of the spatially explicit three-patch solutions. We derive the appropriate boundary conditions in Appendix B and show that the leading eigenvalues of the linearized problem cannot be obtained analytically. This approach is therefore not useful in terms of obtaining a spatially implicit approximation to the original PDE problem.

Methods

PDE solutions

The PDE models, Eqs. 1a and 1b with Eqs. 4a–7b, were solved on a three-patch domain with a central good patch as shown in Fig. 1. We used homogeneous Dirichlet boundary conditions and uniform initial conditions. The domain size was kept constant at 40 units, and the good patch size L was varied from 1 to 19 units. We define b i as the patch size of the ith bad patch, and L as the size of the central good patch. The patch orientation was chosen as [b1 L b2], where b1 had size \(\lceil\frac{40-L}{2}\rceil\), and b2 had size \(\lfloor\frac{40-L}{2}\rfloor\). The bad patches bordering the single good patch could not be set to the same size since they needed to be integer valued for the numerical solver. The difference in size between the two bad patches however, was at most 1. Since the smallest bad patch sizes were 10 and 11, the difference between b 1 and b 2 was on the order of 10% of the size of the patch. We could therefore assume that our domain was approximately symmetrical. We produced numerical simulations for a smaller domain size of 20 (results not shown) to ascertain that domain size was not affecting the results.

Fig. 1
figure 1

Plot of a domain with alternating good and bad patches. This is the domain on which the PDE system is solved. The white region denotes a good patch while the black regions denote bad patches. This domain is bordered by lethal boundaries (homogeneous Dirichlet boundary conditions)

Full size image

The matlab solver pdepe was used to obtain numerical solutions. With this solver, the equations are spatially discretized, and then implicit multi-step methods with adaptive time-stepping are used to solve the resulting ODEs (Shampine and Reichelt 1997). Numerical solutions converged to cyclic predator–prey solutions, as expected, and we confirmed that the cycles were indeed periodic in time and that the choice of initial conditions does not affect the results of the spatially explicit RM, May, and VT models. We chose to measure two dynamics of the population, namely average population density and population cycle amplitude. Due to the symmetry of the problem, the population densities are maximized in the centre of the patch and so we measured the population density and cycle amplitude there. Both quantities were measured after a sufficiently long time interval to ensure that transients had disappeared. The calculation of average densities was produced by averaging the maximum and minimum population densities, and the cycle amplitude was obtained by taking the difference between the maximum and minimum densities.

ODE solutions

The ODEs (Eqs. 18a–21b) were solved using the Matlab ODE solver ODE45, which is a pair of Dormand-Price 4th- and 5th-order solvers (Shampine and Reichelt 1997). The ODE models were tested for good patch sizes 1 ≤ L ≤ 20. Average population densities and cycle amplitude were calculated using the same method as given in the PDE solution methods.

Comparison of ODE and PDE results

Our simulation results are shown in Fig. 2 (cycle amplitude) and Fig. 3 (population average densities). Figure 4 shows the relative error between the spatially explicit and implicit models. These figures illustrate that habitat fragmentation produces similar trends in the solution behaviour of the spatially explicit and implicit models. Specifically, population dynamics under the ODE and PDE frameworks have similar trends in population average for larger good patch sizes, L. For L ≤ 7 however, the two models behave somewhat differently. The relative error between the spatially implicit and explicit models is less than 25% for average densities and less than 35% for cycle amplitude for L > 7. This is with the exception of the relative error in cycle amplitude when using LV reaction terms, which will be discussed later.

Fig. 2
figure 2

Plots of amplitude (scaled between zero and one) of predator and prey against good patch size using the LV, RM, May, and VT reaction terms with habitat fragmentation in the ODE and PDE models. The PDE figures are for a domain with patch sizes [b 1 L b 2], where L = 1 − 19, \(b_1=\lceil\frac{40-L}{2}\rceil\), and \(b_2=\lfloor\frac{40-L}{2}\rfloor\). The figure shows prey and predator amplitude, for the ODE (a, c) and PDE (b, d) models

Full size image
Fig. 3
figure 3

Plots of average (scaled between zero and one) of predator and prey against good patch size using the LV, RM, May, and VT reaction terms with habitat fragmentation in the ODE and PDE models. The PDE figures are for a domain with patch sizes [b 1 L b 2], where L = 1 − 19, \(b_1=\lceil\frac{40-L}{2}\rceil\), and \(b_2=\lfloor\frac{40-L}{2}\rfloor\). The figure shows prey and predator amplitude, for the ODE (a, c) and PDE (b, d) models

Full size image
Fig. 4
figure 4

Plots of relative error in average and amplitude for predator and prey in the ODE and PDE models. The relative error is plotted against good patch size using the LV, RM, May, and VT reaction terms with habitat fragmentation in the ODE and the PDE models. The PDE figures are for a domain with patch sizes [b 1 L b 2], where L = 1 − 19, \(b_1=\lceil\frac{40-L}{2}\rceil\), and \(b_2=\lfloor\frac{40-L}{2}\rfloor\). The figure shows the relative error in prey and predator amplitude (b, d) and average (a, c). The relative error is measured as the absolute value of the difference between the spatially implicit and explicit models divided by the maximum value from the two models and multiplied by 100

Full size image

For small patch sizes, the most significant difference in behaviour between the ODE and PDE models is that the ODE models have a larger critical patch size for the prey and predators. The spatially explicit VT and RM models appear to reach a critical patch size for the prey (VT) and predators (VT and RM) at L = 1 (Fig. 3). The other prey and predator populations in the spatially explicit models do not reach a critical patch size for any L ≥ 1. In contrast, the critical patch sizes for prey and predators in the spatially implicit models are reached for L > 1.

The difference between the ODE and PDE models is due to the construction of the domain in the spatially explicit case. The large low-quality patches on either side of the central good patch in the PDE domain do not act as a lethal boundary on the central good patch, which was the boundary condition assumed in our derivation of the ODE model. Since the good patch is not surrounded by a lethal boundary, we expect the critical patch size to be smaller for the PDE models than the ODE models. This is in fact the case, as can be seen in Fig. 3. The ODE models therefore give an upper bound for the critical patch size of the PDE models.

In general, the ODE (for L > 7) and PDE models predict that the average density decreases for predators and prey as fragmentation increases (L decreases). This is true across models with the exception of the models with LV reaction terms, which show an increase in prey average as habitat fragmentation increases.

The changes in cycle amplitude for the ODE and PDE models with RM, May, and VT reaction terms are fairly consistent in that they all decrease as fragmentation increases. A small difference is that the amplitude of the ODE models show sharper decreases than the PDE models, as a consequence of the larger critical patch size in the ODE models. Additionally, the prey cycle amplitude in the spatially implicit May model increases for larger patch sizes while the corresponding amplitude in the spatially explicit model decreases monotonically. The ODE framework also provides an upper bound on the Hopf bifurcation points, which occur for smaller L under the PDE framework.

The spatially explicit and implicit models with LV reaction terms have opposite trends in cycle amplitude. In the PDE model, both prey and predator populations have decreasing cycle amplitude with increasing habitat fragmentation. In contrast, the ODE model predicts an increase in prey and predator cycle amplitude, followed by a sharp decline to zero as L decreases below the prey critical patch size. It is interesting to note that the decrease in cycle amplitude in the PDE model is not monotonic, in fact it shows dramatic fluctuations in cycle amplitude. Furthermore, the amplitude in the LV PDE model does not equal zero for L ≥ 1 and therefore, the LV ODE model, which equals zero for L > 1, provides an upper bound on the Hopf bifurcation points in the spatially explicit model.

Analysis of the ODE models

The benefit of the simplification to ODEs is that we can now use established analytical tools to understand solution behaviour. Below we determine the steady states and the critical patch sizes for the persistence of prey and predator. Bifurcation diagrams are presented for each model with the good patch size, L, as the bifurcation parameter.

Equilibrium values

The stable coexistence steady-state values are shown in Table 2. Steady states were found by solving the ODE models, Eqs. 18a–21b. The spatially homogeneous steady state of the PDE models can be extracted from the ODE models by taking the limit as L tends to infinity.

Table 2 Steady states of the spatially implicit models with LV, RM, VT, and May reaction terms
Full size table

Critical patch size

The inclusion of good patch size in the spatially implicit ODEs gives rise to a critical patch size for prey and predators.

Prey

The critical patch size for prey, L c , is found by determining at what value of L the prey growth rate, r L , becomes zero. At this point, the prey are on the border between persistence and extinction. This occurs when

$$ r_L=0 \Leftrightarrow r-D_n\left(\frac{\pi}{L_c}\right)^2=0 \Leftrightarrow L_c= \sqrt{\frac{D_n}{r}}\pi. $$
(22)

This result can also be obtained by considering the PDE models in a single homogeneous good patch with homogeneous Dirichlet boundary conditions. In this case we assume that f(n,p,x) = f(n,p) in the PDE model (Eqs. 1a and 1b). Our latter assumption is justified by assuming that the growth rate is approximately constant over the single good patch. This ignores the sharp decrease of the growth rate to zero at the edges of the good patch. This decrease however, occurs over a small spatial length. Thus, our assumption may increase the critical patch size by a small amount, but we consider it negligible. Given these assumptions we arrive at the following set of equations

$$ \frac{\partial n}{\partial t} = D_n \frac{\partial^2 n}{\partial x^2} + f(n,p),\label{eq:primary2} $$
(23a)
$$ \frac{\partial p}{\partial t} = D_p \frac{\partial^2 p}{\partial x^2} + g(n,p),\label{eq:primary2-b} $$
(23b)

where the reaction terms, f(n,p) and g(n,p), are shown in Eqs. 4a–7b, with r(x) = r, a constant. In these spatially explicit PDE models (Eqs. 23a and 23b), the critical patch size (Murray 1993; Okubo and Levin 2001) for prey would be,

$$ L_c=\pi \sqrt{\frac{D_n}{f_n(0,0)}}=\pi \sqrt{\frac{D_n}{r}}. $$
(24)

The table of critical patch size values for prey based on the parameters chosen for this study and for all four models are shown in Table 3.

Table 3 Critical Patch Sizes (CPS) and Hopf Bifurcation (HB) Points in the spatially implicit models with LV, RM, May, and VT reaction terms
Full size table

Predators

To find the critical patch size for predators, we must find the domain size at which the predator numerical response becomes positive. This transition occurs when \(\frac{\partial g}{\partial p}(n_0,0)=0\), where n 0 is the equilibrium density of prey in the absence of predators. We assume that the prey are at carrying capacity in the spatially implicit RM, May, and VT models, that is n 0 = kr L . The approximation of n 0 = kr L is approximately correct, although it is an underestimate of the true population size at the centre of a patch for homogeneous Dirichlet boundary conditions, as seen in Fig. 5. Note that there is no predator critical patch size in the spatially implicit LV model since the prey grow exponentially in the absence of predators, and therefore the numerical response will always be positive in the LV equation. Thus, the L c for predators is the same as for prey in the spatially implicit LV model. This could be unrealistic since you would expect the predators, in general, to require a larger habitat to survive than the prey. The critical patch sizes for predators in the spatially implicit RM, May, and VT models are given below. A list of critical patch size values for predators based on the parameters in all models is shown in Table 3. The critical patch size for predator cannot be calculated analytically (see Appendix A), since it involves linearizing around a non-trivial spatially varying prey steady state which could not be expressed as a simple trigonometric function. Critical patch sizes were thus determined numerically.

Fig. 5
figure 5

Plots of prey density against patch size for Eqs. 1 and 2 in the absence of predators with homogeneous Dirichlet boundary conditions for patch sizes L = 2,4,6,8,10. The reaction terms used for the spatially explicit models are RM (a) and May/VT (b). The value of kr L for each patch size L is shown as dotted lines on the graph

Full size image

Rosenzweig–MacArthur/variable territory

It can be shown that the predator numerical response is exactly the same in the spatially implicit RM and VT models. This makes sense since these models are the same except for the additional density-dependent death term in the VT equations, which is zero when we set p = 0. The predator numerical response changes with prey density as

$$ \frac{\partial g}{\partial p}(r_Lk,0) = \frac{\chi c r_Lk}{d+r_Lk}-\delta_L. \label{eq:rmcp} $$
(25)

Setting (25) to zero, we obtain

$$ \frac{\chi c k \big(r-D_n\big(\frac{\pi}{L_c}\big)^2\big)}{d+k\big(r-D_n\big(\frac{\pi}{L_c}\big)^2\big)}=\delta+D_p\left(\frac{\pi}{L_c}\right)^2. \label{eq:rmcp2} $$
(26)

Rearranging Eq. 26, we arrive at the following quartic equation

$$\begin{array}{lll} L^4 (\chi c k r-d\delta-k\delta r)- L^2\pi^2 (D_n(\chi c k-k\delta)\nonumber\\ \quad + D_p(d+kr)) + D_nD_p\pi^4 =0, \end{array}$$
(27)

which can be solved analytically. Considering only the positive solutions, we obtain

$$ L_c = \frac{1}{(2A)^{1/2}}(B \pm (B^2-4AC)^{1/2})^{1/2}. \label{eq:critrm} $$
(28)

where,

(A):

χc k r − dδ − kδ r,

(B):

\(\pi^2 (D_n(\chi c k-k\delta) + D_p(d+kr))\),

(C):

\(D_nD_p\pi^4\).

Using the parameter values for the spatially implicit RM models given in Table 1, we find that the smaller of the two solutions (Eq. 28) gives a value which is smaller than the prey critical patch size. Therefore, only the larger value is relevant. The smaller critical patch size is in fact a positive equilibrium of the predator, but it corresponds to a negative equilibrium for the prey, as shown in the bifurcation diagram in Fig. 6. Therefore, this equilibrium is not biologically relevant.

Fig. 6
figure 6

Bifurcation plots for prey and predator populations in the spatially implicit LV (a, c) and RM (b, d) models with good patch size, L, as the bifurcation parameter. Note that the small positive stable branch of the spatially implicit RM model before the CPS of the predator. This is a real positive branch but it corresponds to a negative stable steady state for the prey, therefore, the predator does not actually persist in this region

Full size image

May

Repeating the previous calculations for the May model we obtain

$$ \frac{\partial g}{\partial p}(n_0,0) = s_L. \label{eq:mcp} $$
(29)

Setting (29) to zero and rearranging, we obtain

$$ s= D_p\left(\frac{\pi}{L_c}\right)^2. \label{eq:mcp2} $$
(30)

Solving (30) for the critical patch size, we find

$$ L_c= \sqrt{\frac{D_p}{s}}\pi. $$
(31)

In contrast to the critical patch size for the spatially implicit RM and VT models, we find that the critical patch size for the spatially implicit May model is independent of k and r. This calculation assumes that n 0 ≠ 0. If the critical patch size for predators is smaller than that of the prey, then the critical patch size is the same for prey and predators, and thus the predator critical patch size depends upon the prey growth rate, r. Having the prey and predator critical patch size the same only occurs if \(\frac{D_p}{s} \leq \frac{D_n}{r}\). This means that the ratio of predator diffusion to its growth rate is smaller than the ratio of prey diffusion to its growth rate. This is unlikely to occur in nature, unless the habitat size required by the predator is indeed smaller than that of the prey.

Single-patch Robin boundary conditions

We mention here that we can numerically determine the critical patch sizes for prey in the single-patch PDE model with Robin boundary conditions, derived in Appendix B. The details are shown in Appendix B.3. We find that the single-patch Robin boundary conditions approximation yields critical patch sizes that are an underestimate, rather than an overestimate, of the actual critical patch sizes for the three-patch model.

Bifurcation diagrams

We investigated bifurcations in the solution behaviour as the good patch size, L, is varied. This allowed us to gain a better understanding of the steady-state stability and average densities. It also confirmed that the onset of oscillations was due to a Hopf bifurcation. Additionally, the bifurcation plots gave us more precise numerical values of the critical patch size and the Hopf bifurcation points since the ODE simulation using Matlab assumed the patch size was integer valued. The bifurcation diagrams were created using XPPAUT (Ermentrout 2002). A stiff adaptive solver was used with an error tolerance of 0.0001 and a maximum/minimum allowable step size of 0.1/0.00001. The spatially implicit LV and RM models have bifurcation diagrams shown in Fig. 6 and the spatially implicit May and VT models have bifurcation diagrams shown in Fig. 7.

Fig. 7
figure 7

Bifurcation plots for prey and predator populations in the spatially implicit May (a, c) and VT (b, d) models with good patch size, L, as the bifurcation parameter

Full size image

Since the spatially implicit LV model exhibits oscillations about a centre, it does not give rise to a single Hopf bifurcation but has multiple neutrally stable orbits. In this model, increasing fragmentation (decreasing L) results in an increase in prey density and a decrease in predator density. As fragmentation continues to increase, the prey population eventually reaches extinction at the critical patch size, L c . As discussed in “Critical patch size”, the predator and prey go extinct at the same critical patch size. In contrast to the spatially implict LV model, the spatially implicit RM, VT, and May models all exhibit limit cycle behaviour for large enough L, and the bifurcation diagrams are similar. As habitat fragmentation increases in these models, the cycle amplitude decreases, except for the prey in the spatially implicit May model. This is consistent with the ODE simulations shown in Fig. 2. Upon increasing fragmentation, the system is stabilized by passing through a Hopf bifurcation, after which the predator population decreases until it reaches extinction at the predator critical patch size. As the predator population decreases, the prey population increases until reaching a maximum average density at the predator critical patch size. As L is decreased further, the prey population decreases until it reaches extinction at the prey critical patch size. All of the numerical values of prey and predator critical patch sizes as well as Hopf bifurcation points (Table 3) were verified using MAPLE (results not shown). These values are consistent with the numerical simulations (Figs. 2 and 3).

Nullcline behaviour

In order to understand the underlying model structures behind the observed behaviours, we plotted the prey and predator nullclines for three different patch sizes (Fig. 8). Patch size values for the spatially implicit RM, May and VT models were chosen so that coexistence steady states at the highest and lowest patch sizes correspond to the steady state as an unstable and stable focus, respectively. The middle patch size was chosen so that the steady state was close to the Hopf bifurcation point. As habitat fragmentation increases, the steady state changes from an unstable focus to a stable focus. Thus, habitat fragmentation can stabilize oscillatory dynamics.

Fig. 8
figure 8

Plots of predator and prey nullclines as patch size decreases. The nullclines are labelled with the good patch size, L for the spatially implicit LV (a), RM (b), May (c), and VT (d) models. The nullclines of the prey equation are shown in solid lines while the nullclines of the predator equation are shown in dashed lines

Full size image

The nullcline diagrams show the effect of habitat fragmentation (decreasing L) in the phase plane. The coexistence steady state, which is given by the intersection of the two nullclines, is shown for each model. All models shift to lower predator densities as good patch size decreases. The spatially implicit LV, RM, and May models shift to higher prey densities, while the spatially implicit VT model shifts to lower prey densities. The spatially implicit May model does not consistently shift to lower predator densities but increases in predator density for the middle patch size. These nullcline plots are consistent with the bifurcation diagrams in “Bifurcation diagrams” (Figs. 6 and 7).

The nullclines shift in predictable ways to give rise to these changes in the coexistence steady state. In the spatially implicit LV model, the prey nullcline (solid line) shifts down (predator average increases) and the predator nullcline (dashed line) shifts to the right (prey average increases) as patch size decreases. In the spatially implicit RM, VT, and May models, the prey nullcline shifts down and to the left while keeping the same intercept on the predator axis, as patch size decreases. This behaviour is consistent with the notion that as patch size decreases, the carrying capacity of the prey population also decreases and the patch will thus support smaller populations of prey. In the spatially implicit RM model, the predator nullcline shifts to the right with decreasing patch size. The spatially implicit May and VT models have predator nullclines which are lines passing through the origin with some slope m. As the patch size decreases, the slope m decreases.

Discussion

In their original paper on habitat fragmentation, (Strohm and Tyson 2009) explored the response of a series of predator–prey models, namely Eqs. 1a and 1b with Eqs. 4a–7b, to increased habitat fragmentation. Their approach was chiefly numerical since analysis of the nonlinear PDE models was not possible. The mathematical reasons behind the observed behaviour were therefore not presented, and the authors could only confirm that their results were consistent with published biological investigations. In this paper, we have successfully used spatially implicit ODE models to explain the previously observed effects of habitat fragmentation in the PDE models. Here, we summarize the main points of the paper and discuss the insights obtained through our analysis.

The simplified models The spatially implicit ODE models were derived as simplifications of the original systems in two steps: first, the spatially explicit three-patch system of a good patch flanked by bad patches was approximated with the same equations in a single patch with homogeneous Dirichlet boundary conditions. Second, we linearized the PDE model around the extinction steady state and a spatially uniform approximation to the prey-only steady state and obtained an ODE model based on the principal eigenvalues of the two states. The resulting ODE models are thus only roughly related to the original PDE models, and yet they are very useful tools in the interpretation of the PDE solutions.

Our results were obtained using an alternate form for the logistic growth term in Eqs. 4a–7b. The choice of the logistic function, as found in Eq. 3, with the growth rate inside the logistic term, resulted in a functional dependence between the patch size, L, and both the growth rate and carrying capacity of prey (and predator in the spatially implicit May model). Biologically, this means that smaller patch sizes can support fewer individuals. In contrast, if the standard form of the logistic equation (Eq. 2) is used, a change in the growth rate does not also change the carrying capacity. This behaviour is problematic since it implies that the carrying capacity of each population is independent of patch size. The alternate logistic term thus makes more biological sense in the context of the current work. Nevertheless, since the standard form of the logistic term is historically used, we investigated the effect of changing the growth term to the standard form (unpublished work) and, while there were some differences in the results, the overall trends were very similar.

Comparison of the ODE and PDE solution behaviour We emphasize here that the ODE models are obtained from two extreme simplifications of the system—three patches to one patch with Dirichlet boundary conditions and linearization around the extinction steady state. One would therefore not expect the simplified system to have much bearing on the solution behaviour of the original three-patch system for parameter values where the populations coexist. We have shown, however, that the ODE models are nonetheless very useful in terms of developing our understanding of the coexistence solution behaviour of the spatially explicit system. The ODE models were able to reproduce many of the same trends as the PDE models as patch size decreases: initially, decreases in L do not affect the population averages or cycle amplitudes. As L decreases further, there is a precipitous drop in cycle amplitudes and in predator densities at approximately the same value of L. In two of the PDE models, as in the ODE models, the prey density decreases slightly then increases as predator populations drop, and then finally drops precipitously toward zero as L decreases.

For larger good patch sizes L, the match between the PDE and ODE models is very good. Differences in the models appear as patch size decreases and begins to have a stronger effect on the solutions. In particular, the predator densities drop more quickly in the ODE models than in the PDE models. Consequently, the ODE prey densities increase as L decreases, as a result of the rapidly decreasing predator populations. The prey densities in the LV and RM PDE models show the same behaviour, though the increase in prey density is much less marked owing to the more gradual decrease of predator densities. The prey densities in the May and VT models simply decrease.

In contrast to the spatially implicit models, nearly all of the spatially explicit models did not reach a critical patch size L c within the patch sizes L tested in our study. If a critical patch size was reached, it was only for the lowest value of L studied, and thus the extinction predictions of the spatially implicit models do not match those of the spatially explicit models for small L. This is consistent with the study of Cobbold et al. 2005 who found that their derivation of a reduced spatially implicit model only holds under the assumption that the patch sizes are sufficiently large. Since the critical patch sizes we obtained are consistently larger for our ODE models than for our PDE models, we conclude that the ODE models provide an upper bound on the critical patch size of the PDE models. We would expect that the critical patch size in the PDE model would increase as the rate of diffusion of the species increases. This is due to the construction of the domain, since higher diffusion rates would cause a greater density of prey and predator to be lost at the Dirichlet boundaries. Therefore, we would expect the results of the ODE and PDE match more closely as diffusion is increased.

By inspecting the bifurcation diagrams (Figs. 6 and 7) with respect to patch size for the RM, VT, and May models, we observe a series of transitions in solution behaviour. As fragmentation increases, the system exhibits cyclic predator–prey dynamics, then stable predator–prey populations, then extinction of the predator, and finally, extinction of the prey. In the case of the LV model, predator and prey extinction occurs simultaneously. In terms of cycle amplitude, the RM, May, and VT ODE and PDE models responded similarly to habitat fragmentation. In all models studied, the onset of cycles due to a Hopf bifurcation occurred for smaller patch sizes in the PDE models than in the ODE models. Thus, the ODE models provide an upper bound for both the prey and predator critical patch size and the Hopf bifurcation points. The precipitous drop in cycle amplitude and population density is also apparent in the bifurcation diagrams.

The spatially implicit and explicit models differ significantly in the calculation of amplitude using LV reaction terms (Fig. 4). This is likely due to nonlinear effects of the patch size, L, on the LV reaction terms. Initial conditions and nonlinear effects are known to change the neutrally stable orbit that the LV model solutions will follow. Therefore, caution should be exercised if using this technique on models with neutrally stable orbits around a centre.

New mechanistic insights The analysis of the ODE models provides some interesting insights into the biological system under the effects of habitat fragmentation. At low levels of fragmentation, the system still allows persistence of both prey and predator but with some decreased abundance. If the patch size is sufficiently reduced, the dynamics of the predator–prey cycle will stabilize. Around the critical patch size, there is a sudden decrease in average population density. Therefore, based on the ODE results, the effects of increasing fragmentation may not be realized until the populations are at a very fine balance between persistence or extinction. From a management standpoint, this behaviour can pose problems since the decrease in predator and prey abundance is relatively slow up to the critical patch size for predators and prey.

The overall trends in average densities for all models are consistent with the results of Fahrig 2003 and Ryall and Fahrig 2006. Some studies they reviewed found that, similar to the RM, VT, and May models, decreasing good patch size had consistently negative effects on biodiversity and the abundance of predators and their prey. Other studies found that, similar to the LV models, decreasing good patch size causes decreases in the abundance of predators, but can result in an increase in the abundance of their prey.

The loss of oscillations in population densities due to habitat fragmentation in the spatially implicit and explicit models shows that habitat fragmentation stabilizes the dynamics of the system. Stabilization of predator–prey systems can also occur if disease is included in the predator population (Hilker and Schmitz 2008). Predator disease increases the death rate of predators, thereby making them more dependent upon the prey to survive (Hilker and Schmitz 2008). In our study, since the death (birth) rate of the predators increases (decreases) with increasing fragmentation, a similar increasing dependence upon prey may be stabilizing the population dynamics of the system.

Typically, diffusion rates for predators are larger than those for prey, i.e. D p  > D n . With the modified growth and death rates, Eqs. 16a and 16b, we therefore arrive at the result that predator death rates will increase more quickly than prey growth rates will decrease as patch size L decreases. Relating this result to the spatially explicit situation of predator–prey interactions in fragmented habitat, we find that the result is consistent with the notion that predators disperse further than prey over the same time interval. Thus, if patch size decreases, predators are more likely to disperse past the edges of a good patch and into neighbouring bad patches. Therefore, they should be more strongly affected by habitat fragmentation. This behaviour is confirmed by our analysis, and is supported by the studies of Ryall and Fahrig 2006. They found that habitat fragmentation had greater negative effects on the abundance of specialist predators than on their prey. In another study by Cobbold et al. 2005, the effect of habitat fragmentation was influenced by the dispersal distance of the species in question.

The ODE model derivation shows a clear relationship between habitat fragmentation and population growth rate. Our work suggests that fragmentation of habitat acts by effectively decreasing the growth rate of the prey and increasing the death rate of the predator (in the spatially implicit May model, it is the predator growth rate rather than death rate that is affected by habitat fragmentation). These functional relationships are consistent with the idea that in fragmented habitats it may be difficult for prey to reproduce due to lack of resources. Furthermore, decreased densities of prey will result in decreased predator densities.

We can develop additional insights into the biological system from the nullcline plots of the ODE models in Fig. 8, as these reveal the underlying structure of the dynamics. In particular, the movement of the coexistence steady-state point as a function of changes in the growth and death rates (Eqs. 16a and 16b), as good patch size L decreases, serves to elucidate why the predator–prey solutions transition from cycles to steady state and, ultimately, extinction as fragmentation is increased. The RM, May, and VT models all have very similar nullcline structure, and the coexistence steady state moves from one side to the other of the prey nullcline maximum as good patch size L decreases. Fundamentally therefore, the response of the different models to habitat fragmentation is governed by the same underlying structure. These three models therefore, have equivalent value in terms of providing explanations for the response of the PDE models to habitat fragmentation. From a biological perspective, the VT model is the most realistic of the three, but the additional complexity does not affect the basic dynamics. This result is an argument in favour of reliance upon the simpler RM model which is more analytically tractable. The LM model is the simplest of the four, but it is clearly governed by a very different underlying structure, and so the simplifications inherent in the LM equations are too drastic to yield reliable predictions of the PDE solution behaviours.

Further considerations The main differences between the spatially explicit and implicit models are due to the construction of the domain in the PDE models. The PDE models incorporated spatial heterogeneity through a spatially varying growth rate, which was considered zero in low-quality patches. Other authors have assumed that the growth rate could be negative in low-quality patches (Cantrell and Cosner 1991). In this case, the low-quality patches surrounding a central high-quality patch would act more like a lethal boundary than in the current model framework. Assumption of the lethal boundaries was key in the simplification to the ODE model, and therefore, we would expect the ODE and PDE models to produce more similar results if low-quality patches were characterized by negative growth rates.

Our use of zero rather than negative growth rates in bad patches in the original PDE models was motivated by the observation that populations can generally disperse into low-quality patches without experiencing instant mortality (lethal boundaries), for instance, in the case of a forest clearcut. The boundaries to a high-quality patch are realistically not reflecting boundaries either, unless we are dealing with an island or a habitat surrounded by a fence, such as a nature reserve. Therefore, we would expect the good patch to be bordered by boundaries that are neither entirely reflecting or absorbing, but most likely some combination of the two. This observation suggests that simplifying the PDE models using partially absorbing (Robin) rather than Dirichlet boundary conditions may result in a spatially implicit ODE model with solution behaviour that is a closer match to the PDE solution behaviour. This reasoning is consistent with the results of our study. We found that the Robin boundary condition approximation resulted in population average and population cycle amplitude dynamics which more closely followed those of the multi-patch PDE system. The Robin boundary condition approximation provided a lower bound on the critical patch size and Hopf bifurcation points.

Future work Our work has approximated the multi-patch PDE model with a single-patch PDE model, then reduced it to the appropriate spatially implicit model. Future work could be done to extend this analysis by deriving the principal eigenvalue for the multi-patch PDE system using the theory of Fourier series while including the spatially varying prey growth rate. An additional level of complexity would be to include the spatially varying advection, which was included in the original model (Strohm and Tyson 2009). The resulting problem would not be self-adjoint, so the reduction to ODEs could not be done in the same way, using Fourier series, but the principle eigenvalue of the spatial operator (Keener 2000; Kielhofer 2004) would still determine the bifurcation structure of the PDE.