The effect of habitat fragmentation on cyclic population dynamics: a reduction to ordinary differential equations

Abstract

Habitat fragmentation is known to be a key factor affecting population dynamics. In a previous study by Strohm and Tyson (Bull Math Biol 71:1323–1348, 2009), the effect of habitat fragmentation on cyclic population dynamics was studied using spatially explicit predator–prey models with four different sets of reaction terms. The difficulty with spatially explicit models is that often analytical tractability is lost and the mechanisms behind the behaviour of the models are difficult to analyse. In this study, we employ a simplification procedure based on a Fourier series first-term truncation of the spatially explicit models Strohm and Tyson (Bull Math Biol 71:1323–1348, 2009) to obtain spatially implicit models. These simpler models capture the main features of the spatially explicit models and can be used to explain the dynamics observed by Strohm and Tyson. We find that the spatially implicit models and the spatially explicit models produce similar responses to habitat fragmentation for larger high-quality patch sizes. Additionally, we find that the critical patch size of the spatially implicit models provides an upper bound on the critical patch size of the spatially explicit models. Finally, we derive an approximation of the multi-patch habitat by a single-patch habitat with partial flux boundary conditions which allows for a lower bound on the critical patch size to be calculated.

Appendices

Appendix A: Linearization around the prey-only steady state

A.1 Linearization about prey-only steady state

Instead of assuming that prey have a fixed, uniform population size, we investigate here the linearization of the model around the true spatially varying prey-only steady state. So, consider the full PDE model (Eqs. 1a and 1b) with Rosenzweig–MacArthur reaction terms (Eqs. 6a and 6b) and homogeneous Dirichlet boundary conditions.

$$ \frac{\partial n}{\partial t} = D_n \frac{\partial^2 n}{\partial x^2} + \left( r-\frac{n}{k}\right)n -\frac{cnp}{d+n},\label{eq:rmpde-a} $$

(32a)

$$ \frac{\partial p}{\partial t} = D_p \frac{\partial^2 p}{\partial x^2} + \frac{\chi cnp}{d+n}-\delta p.\label{eq:rmpde-b} $$

(32b)

Linearizing this set of equations about the (n*,p*) = (0,0) steady state, with n = n* + εn ₁ and p = p* + εp ₁ (ε < < 1) we have the O(ε) equations:

$$ \frac{\partial n_1}{\partial t} = D_n \frac{\partial^2 n_1}{\partial x^2} + rn_1, $$

(33a)

$$ \frac{\partial p_1}{\partial t} = D_p \frac{\partial^2 p_1}{\partial x^2} -\delta p_1. \label{eq:rmpde1} $$

(33b)

From this equation, we see that the predator equation only has the steady state, p ₁ = 0. The prey equation was solved in “Model simplification”

$$ n_1(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), $$

(34)

Considering only the dominant eigenvalue, we have the prey solution n _1,d :

$$ n_{1,d}(x,t)= A_1 e^{[r-D_n(\pi/L)^2]t}\sin\left(\frac{\pi x}{L}\right), $$

(35)

In the absence of predators (p = 0), allow the prey to grow according to Eqs. 32a and 32b. Note that this prey equation will reach some steady state, $\bar{n}(x)$, since the prey population cannot exceed the carrying capacity rk. Note that due to the nonlinear terms in Eqs. 32a and 32b, which come into effect as the prey population grows away from the extinction steady state, the steady state $\bar{n}(x)$ is not a scalar multiple of n _1,d (Kot 2001). We then linearize these equations around the steady state $(n*,p*)=(\bar{n}(x),0)$ with n = n* + εn ₁ and p = p* + εp ₁ (ε < < 1) we have the O(ε) equations:

$$ \frac{\partial n_1}{\partial t} = D_n \left(\frac{\partial^2 n_1}{\partial x^2}+\frac{\partial^2 \bar{n}}{\partial x^2}\right) \!+\! \left( r-\frac{\bar{n}}{k}\right)n_1 - \frac{\bar{n}n_1}{k}-\frac{c\bar{n}p_1}{d+\bar{n}}, \label{eq:preyeps} $$

(36a)

$$ \frac{\partial p_1}{\partial t} = D_p \frac{\partial^2 p_1}{\partial x^2} + \frac{\chi c\bar{n}p_1}{d+\bar{n}}-\delta p_1. \label{eq:rmpde2} $$

(36b)

Since the predator equation (Eq. 36b) is decoupled from the prey equation (Eq. 36a), we can attempt to solve for the eigenvalues of Eq. 36b. Following the results of Cantrell and Cosner 2003, we transform (36b) into the eigenvalue problem

$$ D_p \nabla^2 \psi + m(x) \psi = \sigma \psi, \text{ in } \Omega $$

(37a)

$$ \psi = 0,\text{ at } \partial \Omega \label{eq:eigpred} $$

(37b)

where $m(x)=\frac{\chi c\bar{n}}{d+\bar{n}}-\delta$, ψ is the eigenvalue with corresponding eigenfunction σ, and Ω = [0,L]. This problem can then be treated with variational methods to determine the principal eigenvalue, σ ₁. The positivity of the principal eigenvalue determines whether or not the predator can grow or will become extinct (Cantrell and Cosner 2003). By the results of Cantrell and Cosner, we have

$$ \sigma_1 \!=\! \max\limits_{\psi \in W_0^{1,2}(\Omega),\psi \neq 0} \left( \frac{-D_p\int_{\Omega}\mid \nabla \psi \mid^2 \mathrm{d}x + \int_{\Omega}m(x)\psi^2 \mathrm{d}x } {\int_{\Omega}\psi^2 \mathrm{d}x} \right) $$

(38)

where ψ is assumed to be zero at the boundaries of the domain, x = 0,L, and $\mid \nabla W_0^{1,2} \mid^2$ must be integrable over Ω. This equation is very difficult to solve analytically due to the spatially varying prey population at equilibrium, $\bar{n}(x)$, in the absence of predator. For this reason, we cannot solve this equation analytically and must look to the ODE simplification approach obtained by linearization around the extinction steady state and an assumed uniform prey-only steady state to answer questions regarding persistence of predators.

Appendix B: Reduction of the multi-patch system to a single-patch system

B.1 Reduction to a single-patch with partial flux boundary conditions

Consider a three-patch domain (0,l ₃) in which the outer two patches ((0,l ₁) and (l ₂,l ₃)) are “bad” patches and the central patch is a “good” patch. We have homogeneous Dirichlet boundary conditions at x = 0,l ₃. Consistent with the framework of our spatially explicit PDE model, the bad patches are areas in which there is no prey growth. As before (“Model simplification”), we wish to approximate the PDE solutions on this three-patch domain with solutions to an ODE model on a single good patch. This time however, instead of assuming homogeneous Dirichlet boundary conditions at either end of the good patch, we will use the PDE model to derive appropriate boundary conditions for the single-patch ODE model. For simplicity, we base our analysis on the LV model (Eq. 8). Our equations become:

$$ \frac{\partial n_m}{\partial t} = rn_m+D_n\frac{\partial^2 n_m}{\partial x^2}, \quad l_1 \leq x \leq l_2, $$

(39a)

$$ \frac{\partial n_{lr}}{\partial t} = D_n\frac{\partial^2 n_{lr}}{\partial x^2}, \quad 0 \leq x \leq l_1, \quad l_2 \leq x \leq l_3, \label{eq:pfsimple} $$

(39b)

where n _m is the prey population inside the middle good patch and n _lr is the prey population in the bad patches on the left and right of the central good patch. Matching fluxes at the interfaces in between patches and employing boundary conditions, we have:

$$ n_m(l_1,t)=n_{lr}(l_1,t)\quad \text{and} \quad n_m(l_2,t)=n_{lr}(l_2,t),\label{eq:pfbcs-a} $$

(40a)

$$ \frac{\partial{n_m}}{\partial{x}}(l_1,t)=\frac{\partial{n_{lr}}}{\partial{x}}(l_1,t)\quad \text{and} \quad \frac{\partial{n_m}}{\partial{x}}(l_2,t)=\frac{\partial{n_{lr}}}{\partial{x}}(l_2,t),\label{eq:pfbcs-b} $$

(40b)

$$ n_{lr}(0,t)=0 \quad \text{and} \quad n_{lr}(l_3,t)=0.\label{eq:pfbcs-c} $$

(40c)

We take the approach of (Ludwig et al. 1979) and look for steady-state solutions to Eq. 39b that satisfy the conditions in Eqs. 40a–40c. We define

$$ n_{lr}(x,t) = \left\{ \begin{array}{lr} n_l(x,t) & \;\;\mbox{if\;\; $0 \leq x \leq l_1$}, \\ n_r(x,t) & \;\;\mbox{if\;\; $l_2 \leq x \leq l_3$}. \end{array} \right. $$

(41)

Steady-state solutions in the outside bad patches can be easily obtained:

$$ n_{l}(x)=c_1x+c_2,\quad n_{r}(x)=c_3x+c_4, $$

(42)

where l/r denotes the left and right bad patches, respectively. To satisfy the boundary conditions (Eqs. 40a–40c) at x = 0,l ₃, we must have that c ₂ = 0 and c ₄ = − c ₃ l ₃. Taking spatial derivatives of our solutions, we obtain,

$$ \frac{{\rm d}n_{l}}{{\rm d}x}(x)=c_1,\quad \frac{{\rm d}n_{r}}{{\rm d}x}(x)=c_3, $$

(43)

Matching the population density and fluxes at x = l ₁ and x = l ₂ results in:

$$ n_m(l_1)=n_{l}(l_1)=c_1l_1=\frac{{\rm d}n_{l}(l_1)}{{\rm d}x}(l_1)=\frac{{\rm d}n_m(l_1)}{{\rm d}x}(l_1),\label{eq:pfbcs2-a} $$

(44a)

$$\begin{array}{rll} n_m(l_2)&=&n_{r}(l_2)=c_3l_2-c_3l_3\nonumber\\ &=&\frac{{\rm d}n_{r}(l_2)}{{\rm d}x}(l_2-l_3)=\frac{{\rm d}n_m(l_2)}{{\rm d}x}(l_2-l_3),\label{eq:pfbcs2-b} \end{array}$$

(44b)

The conditions (44a, 44b) can be used to set up the boundary conditions between the middle and left/right patches in a number of ways. In particular, it is possible to write them entirely in terms of n _m . We can thus reduce our three-patch problem to a single-patch problem for n _m (x). Our problem becomes (with n _m (x,t) replaced by n(x,t)):

$$ \frac{\partial n}{\partial t} = rn+D_n\frac{\partial^2 n}{\partial x^2},\label{eq:pfsimple2-a} $$

(45a)

$$ \frac{\partial n}{\partial x}(l_1,t) = \frac{n(l_1,t)}{l_1},\label{eq:pfsimple2-b} $$

(45b)

$$ \frac{\partial n}{\partial x}(l_2,t) = \frac{n(l_2,t)}{l_2-l_3}.\label{eq:pfsimple2-c} $$

(45c)

The partially absorbing boundary conditions above (often referred to as homogeneous Robin boundary conditions) have population size scaled by the size of the surrounding bad patches. Note that the boundary condition at x = l ₂ has an additional negative sign, which allows the boundary conditions to be symmetric (assuming that the surrounding bad patches are the same size, l ₁ = l ₃ − l ₂). As the bad patches become smaller, the boundary conditions approach homogeneous Dirichlet boundary conditions and as they get larger the system approaches homogeneous Neumann boundary conditions.

In the three-patch model, Fig. 1, we have fixed the domain size and have positioned the good patch at the centre of the domain. Therefore, the bad patches surrounding the central good patch will have identical size defined by B = l ₁ = l ₃ − l ₂. For convenience, we shift our single good patch in Eqs. 45a–45c by l ₁, so that our domain is defined over (0,L), where l = l ₂ − l ₁. This results in the problem:

$$\frac{\partial n}{\partial t} = rn+D_n\frac{\partial^2 n}{\partial x^2},$$

(46a)

$$ \frac{\partial n}{\partial x}(0,t) = \frac{n(0,t)}{B}, \label{eq:bc1} $$

(46b)

$$ \frac{\partial n}{\partial x}(L,t) = \frac{n(L,t)}{-B}, \label{eq:bc2} $$

(46c)

Solving this system of equations using the Fourier series approach as before, we find that the dominant eigenvalue for prey growth is given by

$$ \tan(\gamma_m L) = \frac{2\gamma_m B}{\gamma_m^2B^2-1} \label{eq:trans} $$

(47)

where $\gamma_m =\sqrt{\frac{-\mu}{D_n}}$ for m = 0,1,.... Solving this equation for μ, we then obtain the dominant eigenvalue

$$ \lambda_R=r+\mu. \label{eq:RobinDomEig} $$

(48)

Since Eq. 47 cannot be solved analytically, we cannot obtain a closed form expression for the leading eigenvalue, λ _R . It is therefore not instructive to attempt to simplify to an ODE based on this leading eigenvalue as we did for the single-patch case with homogeneous Dirichlet boundary conditions.

Even though (47) can not be solved analytically, it can be solved numerically. Figure 9 shows that as the patch size, L, decreases, the dominant eigenvalue decreases slowly at first, and then rapidly for small patch sizes. Changes in good patch size thus have the strongest effect for small values of L. This is consistent with the results we obtained for a single patch with Dirichlet boundary conditions, which is also shown in Fig. 9. We also notice that the dominant eigenvalue is larger when the boundary conditions are Robin rather than Dirichlet. We point out that the Robin boundary condition eigenvalue is strictly positive for the patch sizes tested numerically. In Appendix B.3, we compute the critical patch size for the Robin boundary condition case, and so the dominant eigenvalue for each model does eventually become negative, but at a much smaller value of good patch size than for the Dirichlet boundary condition case.

Fig. 9

Plot of the dominant eigenvalue as a function of good patch size, L for the single-patch PDE system with Dirichlet (a) and Robin (b) boundary conditions using LV, RM, May, and VT reaction terms. In this figure, we vary L from 1 to 20 and compute bad patch size B using $B=\frac{40-L}{2}$

B.2 Comparison of PDE results for multi-patch and single-patch systems

In the paper, we compared the solutions of the multi-patch model with homogeneous Dirichlet boundary conditions to the solutions of the ODEs derived from a single-patch approximation with Dirichlet boundary conditions on the single patch. Here, we compare the multi-patch solutions to single-patch solutions when the single patch has Robin boundary conditions. These are defined in Eqs. 46b and 46c. Since an analytic solution to Eq. 47 is not possible, we cannot analytically reduce the single-patch PDE system to an ODE, but we can compare the PDE solutions obtained on the two different types of patches (multi-patch with Dirichlet BCs and single-patch with Robin BCs).

Our simulation results for the multi-patch PDE model with homogeneous Dirichlet boundary conditions and the single-patch PDE model with homogeneous Robin boundary conditions are shown in Figs. 10 and 11. The single-patch PDEs were solved using the same matlab solver and methods as the multi-patch PDE models in “PDE solutions”, with the exception that the domain was a single good patch of length L varied from 1 to 19 units, and the boundary conditions were of the homogeneous Robin type. The value of B, the bad patch size, in the specification of Eqs. 46b and 46c varies with L and was defined by $\frac{40-L}{2}$. The single-patch PDE model with homogeneous Robin boundary conditions has similar trends in amplitude and average as the multi-patch PDE model, with a few exceptions. The general trends are that the amplitude and average decrease as the habitat fragmentation increases. The single-patch and multi-patch PDE models differ since the averages and amplitudes in the single-patch PDE model decrease more slowly than in the multi-patch model. Additionally, the average in the single-patch PDE model with RM reaction terms does not have a small increase around L = 3,4 which occurs in the multi-patch PDE model. In terms of amplitude, the major difference is that the multi-patch spatially explicit PDE models with May, RM, and VT reaction terms have a Hopf bifurcation for L ≥ 1, whereas the single-patch PDE models with the May, RM, and VT reaction terms do not. A final difference between the single-patch and multi-patch PDE models is that the amplitude, when using LV reaction terms, show opposite trends. In the multi-patch model, the amplitude decreases while in the single-patch model the amplitude increases. Since the single-patch PDE model exhibits a slower decrease for both amplitude and average (than the three-patch PDE model), we conclude that this approximation may provide a lower bound on the critical patch size of the multi-patch PDE model.

Fig. 10

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 PDE models with partial flux (pf) and homogeneous Dirichlet boundary conditions. The PDE figures with homogeneous Dirichlet boundary conditions are shown for a domain with patch sizes [b ₁ L b ₂], where L = 1 − 19, $b_1=\lceil\frac{40-L}{2}\rceil$, and $b_2=\lfloor\frac{40-L}{2}\rfloor$. The PDE figures in a single patch with partial flux boundary conditions are shown for domain sizes ranging from L = 1 − 19. The figure shows prey and predator amplitude, for the single-patch PDE model with partial flux boundary conditions (a, c) and the three-patch PDE model with homogeneous Dirichlet boundary conditions (b, d)

Fig. 11

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 PDE models with partial flux (pf) and homogeneous Dirichlet boundary conditions. The PDE figures with homogeneous Dirichlet boundary conditions are shown for a domain with patch sizes [b ₁ L b ₂], where L = 1 − 19, $b_1=\lceil\frac{40-L}{2}\rceil$, and $b_2=\lfloor\frac{40-L}{2}\rfloor$. The PDE figures in a single patch with partial flux boundary conditions are shown for domain sizes ranging from L = 1 − 19. The figure shows prey and predator average, for the single-patch PDE model with partial flux boundary conditions (a, c) and the three-patch PDE model with homogeneous Dirichlet boundary conditions (b, d)

B.3 Analysis of single-patch PDE model with partial flux boundary conditions

Following the approach of Ludwig et al. (Ludwig et al. 1979), a critical patch size for the partial flux boundary conditions can be found by looking at the spatial steady-state solution. Setting the time derivative to zero, and solving for the spatial steady state results in the same transcendental equation as in Eq. 47 but with $\gamma=\sqrt{\frac{r}{D_n}}$, we obtain

$$ \tan(\gamma L) = \frac{2\gamma B}{\gamma^2B^2-1} $$

(49)

Since B is a function of L, we can solve this equation numerically using Matlab. The smallest positive L for which this equation has a solution gives the critical patch size for the single-patch PDE model with partial flux boundary conditions. Using the parameter values from Table 1, we find that the critical patch sizes for the PDE model with partial flux boundary conditions using LV, RM, May, and VT reaction terms are 0.2007, 0.0715, 0.0858, and 0.0858, respectively. As reported in Appendix B.2, the partial flux boundary conditions in a single-patch PDE system may give a lower on the critical patch size of the multi-patch PDE system. If this is the case, then we now have bounds, lower and upper, on the critical patch size for our multi-patch PDE model.