How patch size and refuge availability change interaction strength and population dynamics: a combined individual- and population-based modeling experiment

Overview

To investigate the effects of patch size and habitat complexity (represented by refuge availability) on functional-response parameters, we developed an individual-based allometric predator-prey model (for details, see the supplement for an Overview, Design concepts, Detail protocol, Grimm et al. (2006); Grimm et al. (2010)) to mimic the feeding experiments in the laboratory. We assumed that the maximum feeding rate was driven by mechanical and physiological processes such as chewing and digestion and would not scale with patch size or refuge availability. Therefore, we first investigated the maximum feeding rate without any explicit space properties. Second, we modeled a two-dimensional square area to mimic an explicit patch in which both predator and prey can continuously move. The modeled patch consisted of cells all individuals can enter; however, cells may be marked as refuges preventing predation.

The model processes

The first process applied in the model is prey movement (random walk with randomly chosen direction, 0–2π double precision floating number and allometrically calculated distance). The following processes applied in the model are all decisions and actions of the predator (Fig. 1). First, the digestion of the predator is calculated. Subsequently, the algorithm checks if the predator is handling prey (caught in an earlier time step). If not, and the predator’s gut is full (≥ 60%), it rests (not taking further actions). If the predator is not handling prey and is hungry (gut filling < 60%), the predator moves (random walk, see above). After reaching the new position, the predator investigates if it encounters a prey in the cell. If there is a prey individual in the same cell, it will be attacked. If the attack is successful, another prey item is placed randomly into the grid to keep the prey density constant. The predator starts to handle (chew) prey in the next time step.

Variables and parameters

Most species traits regulating the processes described above follow allometric rules (Kleiber, 1961; Peters, 1983; Brown et al., 2004; Brose, 2010), including velocity, V [cm s⁻¹], of both the predator and the prey (Peters, 1983); and the traits of the predator: gut size, G [mg] (Ibarrola et al., 2012), digestion rate, D [mg s⁻¹] (Ibarrola et al., 2012), handling time, T_h [s] (estimated from Rall et al., 2012, see Supplemental Information) and attack success, S_a [unitless] (Brose et al., 2008; Gergs, 2011):

where v₀, g₀, d₀ and h₀ are constants, a_v, a_g, a_d and a_h are the allometric scalings, and M is the body mass of the corresponding individual. Subscripts, _p and _n indicate predator and prey respectively. We used the widespread generalized Ricker’s function (Persson et al., 1998; Persson & Brönmark, 2002b; Persson & Brönmark, 2002a; Wahlström et al., 2000; Brose et al., 2008; Rall et al., 2011) to describe the scaling of attack success depending on body mass. This function consists of the maximum attack success a₀, predator-prey body-mass ratio, R and its optimum R_opt and a shaping parameter, λ. Predator and prey also possessed some state variables to assist their decision making and activities, i.e., the ‘position’ for all individuals; the ‘gut fullness’ and if the predator is ‘still handling’ and an identifier, ‘prey identity,’ to distinguish between the prey individuals.

Parameters’ range

The cell resolution of the square grid, in which the in silico simulations are conducted is 1 cm × 1 cm. As we intended to mimic laboratory experiments, the walls of the grid are set to ‘wall-boundary condition’ (individuals cannot penetrate the walls). We chose twelve patch sizes ranging from 0.2 m × 0.2 m = 0.04 m² (the size of a standard patch in some terrestrial functional response experiments (Brose et al., 2008; Rall et al., 2010; Vucic-Pestic et al., 2010a; Rall et al., 2011; Vucic-Pestic et al., 2011; Kalinkat, Brose & Rall, 2013) to 100 m² (the size of a field patch (Munyaneza & Obrycki, 1997)). The sizes of each patch were: 0.04 m², 0.16 m², 0.64 m², 1.44 m², 2.56 m², 4 m², 16 m², 36 m², 49 m², 64 m², 81 m², and 100 m². The second independent variable we modeled was prey refuge that served as a surrogate for habitat complexity which preventing feeding. We randomly selected refuge cells on the grid for each simulation in a certain percentage of cells in steps of 5% (5%–75% as the ratio of refuge cells to all cells); see Fig. 2 as a case example. These two independent variables are full-factorially simulated. For each simulation run, the refuge distribution is newly drawn. Those randomly chosen cells do not support any feeding by the predator and therefore act as refuges for the prey. The body masses of the predator and prey were set to 100 mg and 1 mg, a common body-mass ratio for animal predatory interactions, close to the optimal feeding ratio of invertebrates (e.g., Vucic-Pestic et al., 2010b; Rall et al., 2011; Kalinkat et al., 2013). We ran each of the in silico feeding trials for 3,600 steps (representing 1 h). The simulation for estimating the maximum feeding rate was repeated 50 times and each prey density dependent simulation was repeated five times. We simulated prey densities from 2⁰ to 2ⁿ as the density when the predator (only one predator per simulation) is satiated. For example, twenty prey densities from 2⁰ to 2¹⁹ are selected for the patch size of 36 m² and 35% refuge-area ratio. Values for the parameters in allometric equations, Eq. (2), are empirically-based and given in Table 1. These values (Table 1) are derived from the same studies where we derived the formulas. Yet the maximum attack success a₀ is taken as the mean of 5 measurements from Gergs (2011). The optimum predator-prey body-mass ratio is consistent with terrestrial invertebrates from Brose et al. (2008).

Functional response fitting

We first calculated the mean maximum feeding rate for the predator-prey pair. We used a generalized linear model (GLM) assuming that maximum feeding rates follow Poisson distribution as feeding rates were count data of non-negative integers of which the error distribution increases with increasing mean. The statistics were ran in R (R Core Team, 2016), but see chapter 13 in Crawley (2007) for details. Subsequently, we used this mean maximum feeding rate as a fixed parameter in the functional response model (Eq. (1)) to estimate the dependencies of the remaining half saturation density and Hill exponent.

We analyzed the feeding data from IBM models using Real’s functional response, Eq. (1). As there is no well-established scaling relationships of functional-response parameters (half saturation density and Hill exponent) to habitat properties investigated here, i.e., patch size and refuge availability, we preliminarily tested whether the scalings of these functional-response parameters followed a power law or exponential function. To reduce the potential influences of interaction terms (between patch size and refuge availability) which may influence the dependencies of the half saturation density or Hill exponent, we included all interaction terms in the preliminary testing (Zuur et al., 2009). We analyzed in total 16 full models and compared them using the Bayesian Information Criterion (BIC), see Table S3. This analysis revealed that the scalings of half saturation densities with patch size and refuge availability can be best described by a power law and an exponential function, respectively: (3)${\displaystyle N_0=C_{N_0}{A}^{a_{N_0}}{e}^{b_{N_0}R}{e}^{{\gamma }_{N_0}ln\left(A\right)R}}$ where C_N0 is a constant, a_N0 is the scaling exponent of half saturation density to patch size, A, b_N0 is the scaling parameter of half saturation density to refuge availability, R and γ_N0 is the parameter giving the strength of the interaction between patch size and refuge availability. Preliminary analyses also showed that the Hill exponent depended on patch size and refuge availability both following power laws: (4)${\displaystyle h=C_h{A}^{a_h}{R}^{b_h}{e}^{{\gamma }_{h}ln\left(A\right)ln\left(R\right)}}$ where C_h is a constant, a_h is the scaling parameter of the Hill exponent to patch size, A, b_h is the scaling exponent of the Hill exponent to refuge availability, R and γ_h is the parameter giving the strength of the interaction between patch size and refuge availability.

We fitted the functional response model, Eq. (1) with the dependencies described above using a maximum likelihood method, ‘mle2()’ (Bolker & R Development Core Team, 2014), (see Bolker (2008) for details). As we replaced eaten prey after each feeding event (see above), we assumed that the residuals followed a negative binomial distribution. We fitted this functional response model to the data assuming a log link between data and model:

i.e., we did not fit the values for the constants C in Eqs. (3) and (4), but for the intercepts in the ln-transformed version ln(C_N0) and ln(C_h) in Eq. (5). We performed a model selection using the Bayesian Information Criterion (BIC) by comparing all possible combinations of setting the parameters a, b and γ to “0,” resulting in 25 meaningful combinations (note that either a or b only can be excluded if the interaction term, γ is excluded).

Equilibrium densities of the predator-prey system and extinction boundaries

The predator-prey population model has a set of non-trivial analytical solutions, being a predator isocline (8)${\displaystyle N={\left(\frac{x{N_0}^h}{e\omega f_{max}-x}\right)}^{\frac{1}{h}}}$ and a prey isocline (9)${\displaystyle P=r{N}^{1-h}\left(k-N\right)\frac{{N_0}^h+N^h}{k\omega f_{max}}.}$ After obtaining these isoclines, the equilibrium densities of predator and prey are compared with extinction boundaries [# m⁻²]. Such boundary is set to an artificial small number in the common practice of ODE models, but we explicitly set it to two individuals per patch. In cases where the predator population is not sustained, i.e., the equilibrium density is less than the extinction boundary, prey population would grow to its capacity, N = K and the predator population goes extinct, P = 0.

Parameter values for the ODE

The functional-response parameters, the maximum feeding rate, f_max, the half saturation density, N₀ and the Hill exponent, h, are set according to the statistical results of the IBM simulations. We assumed that the predator foraged approximately 12 h a day (Ebeling & Bray, 1976), therefore we added a foraging time proportion $\omega =\frac{1}{2}$. The assimilation efficiency, e, accounts for the proportion of food overwhelmed by the predator which can be converted to its own body mass, which is set to 0.85, a common value for predatory consumers (Yodzis & Innes, 1992; Otto, Rall & Brose, 2007). The prey growth follows the logistic growth consisting of the intrinsic growth rate r and the carrying capacity K. Together with metabolic rate of the predator, these three parameters are calculated by empirically derived equations.

The carrying capacity K scales with body mass, M_n (gram), environmental temperature, T (K), net primary production of the habitat, ${\left({\sigma }_0{e}^{\frac{E_{\sigma }\left(T_0-T\right)}{kT{T}_0}}\right)}^z$, and the trophic level of the prey, tl. The values for all parameters are derived for invertebrate detritivores assuming German weather conditions and productivity: K₀ = e^−31.15; b_K =  − 0.72; E_K = 0.71; k = 8.62e − 05; T = 282.65; σ₀ = 600; E_σ =  − 0.35; T₀ = 293.15; z = 1.03; tl₀ =  − 2.68; tl = 1.5 (see Meehan (2006) and Rall et al. (2010) for details). The growth rate r, scales with body mass (microgram) and environmental temperature, where r₀ = e^32.39, b_r =  − 0.25 and E_r =  − 0.84 (details see Savage et al. (2004a) and Rall et al. (2010)). The metabolic rate x, also scales with body mass (gram) and environmental temperature with x₀ = e^27.68, b_x = 0.72 and E_x = 0.87 (see Peters (1983); Savage et al. (2004b) and Rall et al. (2010) for details). Savage et al. (2004b) reported that field metabolic rate were three times larger than basal, therefore we include the coefficient σ as 3. The normalization constant c_x, $12342.86M_p^{-1}$ (M_p in milligram), converts the metabolism from J s⁻¹ to d⁻¹ (Peters, 1983).

We set predators to 100 mg, and prey to 1 mg, consistent with our individual-based model simulations described above. We also explored the same ranges of patch size and habitat complexity as for the individual-based model simulations explained above. Extinction boundaries for predator and prey were set to two individuals per patch.