Visitors   Views   Downloads
Note that a Preprint of this article also exists, first published June 30, 2016.


The interplay between stability, complexity and biodiversity of ecological networks (especially food webs) is a subject of a long lasting and still ongoing discussion in ecology (e.g., MacArthur, 1955; May, 1972; McCann, 2000). From a mathematical point of view, an increase of species richness is only possible when the interaction strength between the species in a network (i.e., the strength of feeding interactions in a food web) decreases (May, 1972), otherwise species richness (represents biodiversity) must decline. A few possible solutions to overcome this problem have already been proposed, including a non-random organization of the feeding links in real food webs (Yodzis, 1981) and a decrease of specific interaction strength with increasing trophic levels (De Ruiter, Neutel & Moore, 1995), both driven by allometry (Yodzis & Innes, 1992; Brose, Williams & Martinez, 2006; Otto, Rall & Brose, 2007). Classic stability analyses assumed that the strength of interactions (the functional response) increased linearly with increasing resource availability (e.g., May, 1972). In food web studies, interaction strength refers to feeding interactions, which can be studied by investigating the ‘functional response’: how the number of prey eaten by a predator changes with increasing prey densities. In his seminal work, Holling (1959b) described the mechanism of functional responses mathematically and showed that under the simplest assumptions it should follow a hyperbolic curve. This experiment was performed in a simple experimental trial with blindfolded students (the predator) on a plain nine-square foot table as the experimental arena and sandpaper discs as prey. Interestingly, already in the same year Holling (1959a) showed that the shape of functional response could also follow a sigmoid (i.e., s-shaped) curve when investigating small mammals on a large natural scale, including natural complexity in the habitat. The different possible shapes of functional responses are known as type I functional responses (linear with a limit), type II functional responses (hyperbolic) and type III functional responses (sigmoid), forming the core set of functional responses alongside a number of other descendant types (Jeschke, Kopp & Tollrian, 2002). Notably, it has been shown that the type I functional response is an artifact (Jeschke, Kopp & Tollrian, 2004; Sarnelle & Wilson, 2008), so in this study we focus only on the non-linear functional responses. One of the descendants of Holling’s functional response models based on enzyme kinetics (Real, 1977) unifies the type II and the type III functional responses: f N = f m a x N h N 0 h + N h where f(N) is the per capita feeding rate, depending on the resource density, N. The curve is characterized by a maximum feeding rate, fmax [# h−1]; a half saturation density, N0 [# m−2], i.e., the prey density when the predator’s feeding rate reaches half of the maximum feeding rate and a unitless Hill exponent, h, determining the curve shape. If the Hill exponent is set to unity, the resulting function is the “strict” type II functional response. Whereas if the Hill exponent is set to two, it is the “strict” type III functional response in a very classical sense (but see Juliano (2001) for alternative descriptions of the type III functional response). We will subsequently refer to any functional responses as type III functional responses if the Hill exponent is larger than unity.

But why is it so important to know if the functional response is a type II functional response or a type III functional response? The answer is that type III functional responses are known to stabilize population dynamics thereby allowing for persistence of large food webs (Williams & Martinez, 2004; Brose, Williams & Martinez, 2006; Rall, Guill & Brose, 2008), by regulating prey populations to low densities (Nunney, 1980a; Nunney, 1980b). Several mechanisms have been put forward to explain why a type III functional response should appear, including the predator learning to exploit prey better (Holling, 1966) and switching between different prey types to the most abundant prey (Murdoch & Oaten, 1975; Oaten & Murdoch, 1975). More recently, it was suggested that refuges for the prey can also lead to a type III functional response (e.g., Scheffer & De Boer, 1995; Vucic-Pestic et al., 2010a). However, there was only mixed support from studies manipulating habitat complexity in general without introducing prey refuges explicitly (Kaiser, 1983; Hoddle, 2003; Hohberg & Traunspurger, 2005; Hauzy et al., 2010; Vucic-Pestic et al., 2010a; Kalinkat, Brose & Rall, 2013). Those differences may be caused by variations in how habitat complexity influences foraging and hence feeding: (1) complexity negatively affects feeding by e.g., refuges restraining predation especially at low prey densities and eventually leading to a type III functional response; (2) it affects feeding by e.g., obstacles preventing the movements of both predator and prey leading to reduced encounter rates at all prey densities but leaving the functional response type unaffected (Hauzy et al., 2010). Moreover, if the habitat complexity influences neither movement nor refuge provision, it will simply cause a dilution effect, a virtual increase of the patch size (Kalinkat, Brose & Rall, 2013).

Furthermore, functional response studies are predominantly carried out under artificial laboratory conditions (Jeschke, Kopp & Tollrian, 2004; Kalinkat & Rall, 2015). This means that (1) most of the studies mentioned above have used rather artificial habitat complexity and (2) due to spatial limitations of a laboratory, the size of the experimental units are relatively small (e.g., petri-dishes for estimating the functional response of ladybugs, Stethorus japonicus (Gotoh, Nozawa & Yamaguchi, 2004)). Only a few studies have attempted to investigate functional responses in natural environments, but these studies are only roughly comparable to the controlled laboratory studies as they rely on scat counting or gut content analyses combined with assessments of natural prey density  (e.g., Dale, Adams & Bowyer, 1994; Smout & Lindstrøm, 2007). To our knowledge, only one study so far has investigated and compared the simplified laboratory functional response experiments (using petri dishes) with functional responses measured in the greenhouse or in the field (Munyaneza & Obrycki, 1997). In this study, the attack rates in the laboratory were over 40 times higher than those in the greenhouse but those in the field were three to nine times lower than in the petri dishes (note that all functional responses in the original publication were fitted to a type II functional response model and the attack rates were scaled to the total size of the experimental arena, we compared the attack rates scaled to square meters (Rall et al., 2012)). As the experiments carried out by Munyaneza & Obrycki (1997) altered habitat complexity and patch size at the same time, and both gradients increased simultaneously from the petri dish experiments over the greenhouse to the field experiments, it was not possible to disentangle the potential interactive effect of habitat complexity and patch size. Furthermore, Bergström & Englund (2004) reported increases in attack rates with patch size, and studies manipulating habitat complexity reported a decrease in attack rates (e.g., Vucic-Pestic et al., 2010a). This might lead to the explanation that the relatively low attack rates in the field shown in the study of Munyaneza & Obrycki were due to the increased patch size and habitat complexity.

Beside the examples above, we are not aware of any other study addressing the effect of patch size and habitat complexity on the functional-response parameters. Moreover, most of the studies only vary habitat complexity or patch size by up to four levels (e.g., Kalinkat, Brose & Rall, 2013; Bergström & Englund, 2004) and none of them systematically varied both complexity and size. This lack of studies is perhaps due to the laborious nature of functional response studies. For example, Vucic-Pestic et al. (2010b) used prey ranging from one up to 4,000 individuals for fitting a single functional response, measuring up to 90 feeding experiments. Doubling the patch size would already lead to a maximum of 8,000 individuals and a 10 times larger patch would require already a maximum of 40,000 individuals.

As such extreme laboratory settings are not feasible, we developed an individual-based model (IBM) to study the effects of patch size and refuge availability on functional-response parameters. We explored full-factorial patch size and habitat complexity to disentangle effects of both variables and eventually their interactive effects on the functional-response parameters. Subsequently, we analyzed the stability of a predator-prey system depending on patch size and habitat complexity by developing a predator-prey population dynamics model which has a long standing usage and wide applicability (see section “Methods” for details).


Individual-based model of feeding interaction


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.

Schematic diagram of processes in the in silico feeding experiment model.
Figure 1: Schematic diagram of processes in the in silico feeding experiment model.
The text in the dashed oval is the only prey action which is also the start of the processes. Texts in the squares are decisions or actions of the predator. Texts in the diamonds are decisions making, e.g., “handling prey?”

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−1], 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−1] (Ibarrola et al., 2012), handling time, Th [s] (estimated from Rall et al., 2012, see Supplemental Information) and attack success, Sa [unitless] (Brose et al., 2008; Gergs, 2011):

V = v 0 M a v G = g 0 M a g D = d 0 M a d T h = h 0 M p a h , p M n a h , n S a = a 0 R R o p t e 1 R R o p t λ

where v0, g0, d0 and h0 are constants, av, ag, ad and ah 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 a0, predator-prey body-mass ratio, R and its optimum Ropt 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 m2 (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 m2 (the size of a field patch (Munyaneza & Obrycki, 1997)). The sizes of each patch were: 0.04 m2, 0.16 m2, 0.64 m2, 1.44 m2, 2.56 m2, 4 m2, 16 m2, 36 m2, 49 m2, 64 m2, 81 m2, and 100 m2. 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 20 to 2n as the density when the predator (only one predator per simulation) is satiated. For example, twenty prey densities from 20 to 219 are selected for the patch size of 36 m2 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 a0 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).

Schematic diagram of a square grid of the IBM model.
Figure 2: Schematic diagram of a square grid of the IBM model.
The grid in this example is 0.04 m2 with cell resolution of 1 cm × 1 cm and with a ‘wall-boundary condition’ (individuals not able to penetrate the walls). The black cells are the refuge cells. The small white bugs represent prey and the big black bug represents the predator. The random walk of the individuals is decided by a randomly chosen direction, 0–2π and an allometrically decided distance (as denoted by the white circles and black arrows).
Table 1:
Parameters values in allometric equations (Eq. (2)).
Parameter Value Parameter Value
v0 0.546 av 0.29
g0 0.50 ag 0.434
d0 5 × 10−5 ad 0.75
h0 37.504 ah,p −0.330
ah,n 0.173
a0 0.10 Ropt 100
λ 1
DOI: 10.7717/peerj.2993/table-1


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: N 0 = C N 0 A a N 0 e b N 0 R e γ N 0 ln A R where CN0 is a constant, aN0 is the scaling exponent of half saturation density to patch size, A, bN0 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: h = C h A a h R b h e γ h ln A ln R where Ch is a constant, ah is the scaling parameter of the Hill exponent to patch size, A, bh 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:

ln N 0 = ln C N 0 + a N 0 ln A + b N 0 R + γ N 0 ln A R ln h = ln C h + a h ln A + b h ln R + γ h ln A ln R

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(CN0) and ln(Ch) 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).

Population dynamics model

To investigate how patch size, A, and refuge availability (a measurement of habitat complexity), R, affect population dynamics and stability of a predator-prey system, we set up an ordinary differential equations (ODE) model. Such models were widely used to study one population (e.g., Gompertz, 1825; Verhulst, 1838) over food web motifs (e.g., Lotka, 1925; Volterra, 1926; Rosenzweig & MacArthur, 1963; Rosenzweig, 1971; Yodzis & Innes, 1992) to multi-trophic food web models (e.g., Williams & Martinez, 2004; Binzer et al., 2016; Schneider et al., 2016) and used to predict patterns of experimental microcosms (e.g., Schneider, Scheu & Brose, 2012; Fussmann et al., 2014) up to whole food webs (e.g., Boit et al., 2012). The ODE model describes the change in prey density [# m−2], dN, and predator density [# m−2], dP, over time, dt (Rosenzweig & MacArthur, 1963; Yodzis & Innes, 1992; Otto, Rall & Brose, 2007):

d N d t = r N 1 N K ω f m a x N h N 0 h + N h P d P d t = e ω f m a x N h N 0 h + N h P x P

where the prey growth is described by logistic growth with r [d−1] being the intrinsic growth rate and K [# m−2] being the carrying capacity. The prey are consumed by the predator following Real’s functional response, Eq. (1), with fmax being the maximum feeding rate, N0 being the half saturation density and h being the Hill exponent. The predator population grows according to the functional response multiplied by the assimilation efficiency, e, and the effective foraging time proportion ω. Moreover, it loses population density by metabolism, x [d−1].

We used the estimated values from the functional response fitting of our individual-based model (see above) in the ODEs. Additionally we calculated the values for carrying capacity, K, growth rate, r and metabolism, x, according to empirically derived studies (Rall et al., 2010; Meehan, 2006; Savage et al., 2004a; Peters, 1983) (details described afterwards).

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 N = x N 0 h e ω f m a x x 1 h and a prey isocline P = r N 1 h k N N 0 h + N h k ω f m a x . After obtaining these isoclines, the equilibrium densities of predator and prey are compared with extinction boundaries [# m−2]. 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, fmax, the half saturation density, N0 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 ω = 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.

K = K 0 M n b K e E K k T σ 0 e E σ T 0 T k T T 0 z e t l 0 t l 1 r = r 0 M n b r e E r k T x = σ c x x 0 M p b x e E x k T

The carrying capacity K scales with body mass, Mn (gram), environmental temperature, T (K), net primary production of the habitat, σ 0 e E σ T 0 T k T T 0 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: K0 = e−31.15; bK =  − 0.72; EK = 0.71; k = 8.62e − 05; T = 282.65; σ0 = 600; Eσ =  − 0.35; T0 = 293.15; z = 1.03; tl0 =  − 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 r0 = e32.39, br =  − 0.25 and Er =  − 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 x0 = e27.68, bx = 0.72 and Ex = 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 cx, 12342 . 86 M p 1 (Mp in milligram), converts the metabolism from J s−1 to d−1 (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.


Results of IBM simulation

The maximum feeding rate, fmax, was estimated prior to the functional response fitting and revealed that predators of 100 mg fed in average 6.7 prey individuals per hour (Table 2, note that the statistics were performed using a log-link function, i.e., the ln-linear feeding rate was estimated). We subsequently fitted the functional response with a fixed maximum feeding rate. Our model comparison of patch size and refuge availability dependent functional response models, based on BIC, included a scaling of half saturation density with patch size and refuge availability, but it did not include any scaling of the Hill exponent with either parameter (Table 2). The half saturation density increased with refuge availability, and decreased marginally with patch size, see Fig. 3. The estimated Hill exponent across patch sizes and refuge availabilities was 1.284, which is significantly different from a Hill exponent of 1 therefore indicating a type III functional response (Table 2, note that the ln-transformed Hill exponent was tested against “0” what is a Hill exponent of “1”). The emerging functional responses are of the same shape but feeding is realized at higher prey densities with increasing refuge availability (Fig. 4). Other predator-prey body-mass ratios showed similar results, see the section “In silico feeding experiments on other body-mass ratios” in the supplement.

Table 2:
Statistical results for the in silico functional response experiments.
Note that the maximum feeding rate is a priori assumed to be independent of patch size and refuge availability.
Estimate S.E. p-value
Maximum feeding rate ln(fmax) 1.902 0.05 <0.001
Half saturation density ln(CN0) 4.577 0.031 < 0.001
aN0 −0.007 0.005 0.21
bN0 1.777 0.063 < 0.001
γN0 Excluded by model selection
Hill exponent ln(Ch) 0.25 0.011 < 0.001
ah Excluded by model selection
bh Excluded by model selection
γh Excluded by model selection
DOI: 10.7717/peerj.2993/table-2
The effect of patch size (x-axis) and refuge availability (y-axis) on half saturation density (see color scale).

Figure 3: The effect of patch size (x-axis) and refuge availability (y-axis) on half saturation density (see color scale).

Results of the individual-based functional response in silico experiments (black dots) and their corresponding fits (red lines).

Figure 4: Results of the individual-based functional response in silico experiments (black dots) and their corresponding fits (red lines).

The panels are arranged by increasing refuge availability, starting at 5% (A) to 75% (O). The patch size effect on the half saturation density is too small to result in visually distinguishable regression lines. All prey densities are ln-transformed. The grey vertical lines denote the half saturation densities.

Results of population dynamic model

We solved the population dynamics model by a set of analytical solutions (Eqs. (8) and (9)) and the extinction boundaries. In small patches only the prey species survived, but refuge availability relaxed this pattern, allowing predators to survive at smaller patches. Both predator and prey population densities increased with refuge availability, whereas in larger patches, the densities of predator and prey populations decreased slightly (Fig. 5).

Number of surviving species (A), population densities of predator (B) and prey (C) depending on patch size (x-axis) and refuge availability (y-axis).

Figure 5: Number of surviving species (A), population densities of predator (B) and prey (C) depending on patch size (x-axis) and refuge availability (y-axis).

When the system is embedded in very small patches, the predator becomes extinct due to energy limitation and the prey grows to its carrying capacity (green areas in B and C). The non-green areas of B and C show the densities, [# m−2], of predator and prey when the system is feasible (see color code above each plot).


Effects of spatial properties on interaction strength

We developed an allometric individual-based model to investigate the effects of patch size and habitat complexity (represented by refuge availability) on feeding interactions. We found that the interaction strength decreased with refuge availability, as the half saturation density increased with it. This result is consistent with most of empirical studies aiming to account for how refuges affect predation rates (e.g., Kaiser, 1983; Folsom & Collins, 1984; Kalinkat, Brose & Rall, 2013). Our results showed that patch size, however, did not significantly influence the interaction strength. Bergström & Englund (2004) reported that the attack rate (the maximum interaction strength at low prey densities) increased with patch size. This increase was explained by behavioral changes in moving activity (increasing speed of the predator) and an aggregative behavior of both their prey and predator at the walls of their experimental aquariums (animals clustered more at the aquarium walls with increasing patch size). We did not include such behavioral changes in our model as we wanted to provide a simple basic model in this study, which may explain the differences of our results from Bergström & Englund (2004). The Hill exponents in our in silico individual-based experiments were 1.28 across patch sizes and refuge availabilities. This is quite surprising, as a simple type II functional response was thought to be the appropriate model for feeding interaction experiments under simplified conditions in the laboratory. However, the empirical findings of Sarnelle & Wilson (2008) suggested that type III functional responses would emerge if researchers were able to include experimental trials on small prey densities, which was not feasible for experiments carried out in small patches. A few feeding interaction studies on mammals (‘intelligent predators’) carried out in the field also suggested type III functional responses (Holling, 1959a; Smout & Lindstrøm, 2007). More recent studies found type III functional responses for invertebrates as well (Aljetlawi, Sparrevik & Leonardsson, 2004; Vucic-Pestic et al., 2010b). Our study not only corroborates the finding of type III functional response, but also confirms that for the mechanistically simplified predators as in our individual-based model, a type III functional response is appropriate, which is not only suitable for ‘intelligent predators’ with the ability to learn (Holling, 1966). The statistical results for in silico experiments of other body-mass ratios showed consistency with the results discussed above.

Former laboratory experiments that compared a homogeneous habitat with a complex habitat documented a shift from a type II to a type III functional response (Vucic-Pestic et al., 2010a) and argued that this was due to a refuge effect. We did, however, not find any increase in the Hill exponent with increasing refuge availability. As we did not include explicit behaviors for hiding, we infer that this switch from a type II to a type III functional response (or an increase of the Hill exponent) not only needed refuges as shelter for the prey, but also active behavioral changes in sub-habitat choice (Schmitz, Krivan & Ovadia, 2004; Miller, Ament & Schmitz, 2014).

Our individual-based predator-prey model framework allowed us to investigate the effects of patch size and refuge availability on functional-response parameters, which would not have been possible in laboratory or field experiments. Even without incorporating more complex movement models than random walks or behaviors like chasing or hiding, we were able to detect general patterns on the scalings of functional-response parameters with increasing patch size and refuge availability. Nevertheless, future individual based predator-prey models should incorporate more complex movement models to better understand the mechanisms of functional responses.

Effects of spatial properties on population dynamics

To investigate how changes in interaction strength scale up to population dynamics and coexistence, we analyzed a predator-prey ordinary differential equation model. We used the results from our in silico feeding experiments and combined it with empirically measured values for growth, carrying capacity and metabolism (Meehan, 2006; Savage et al., 2004a; Savage et al., 2004b; Brown et al., 2004; Rall et al., 2010). Increasing patch size turned the extinction of predators to survival, meaning the smallest patches were not able to sustain the predator population. This is surprising as the feeding interaction strength does not change with patch size (i.e., a non-significant effect of patch size). As all parameters of the model are constant in respect to patch size, we expected that neither the stability (sensu population dynamics) nor persistence will be affected. This paradox behavior of the system can only be explained by the increasing extinction thresholds with decreasing patch size. We defined the extinction threshold as two individuals per patch leading to increasing densities for extinctions with decreasing patch size (see Fig. 6 as an example). Increasing refuge availability counteracted this pattern and allowed predators to survive at even smaller patches. This is surprising as the half saturation density increased with refuge availability, suggesting less energy intake by the predator. However, increasing half saturation density also led to an increase in prey density in equilibrium that subsequently sustained a higher predator density. Both predator and prey have been feasible in larger patches and across the range of habitat complexities we explored.

Equilibrium densities and extinction boundaries of the predator-prey system.

Figure 6: Equilibrium densities and extinction boundaries of the predator-prey system.

This is an example where 73% of the cells are prey refuges. The solid lines depict the equilibrium population densities of the predator (red) and prey. The dashed line indicates the extinction boundaries of the predator and prey, two individuals per patch.

Using a predator-prey population dynamics model, parameterized by the in silico functional response experiments discussed above, we were able to detect patterns of coexistence when patch size and refuge availability increased. These results are predominantly driven by the incorporation of a nearly constant feeding interactions with increasing patch size and a realistic assumption for the extinction boundaries of populations which is usually ignored in an ordinary differential equation modeling frameworks (e.g., McCann, Rasmussen & Umbanhowar, 2005; Otto, Rall & Brose, 2007). Future studies that aim to investigate effects of space on persistence using ordinary differential equation models should consider to include such more realistic extinction boundaries as presented in our study.


How species interactions react to environmental changes such as habitat homogenization and habitat loss is a key point for understanding how current global changes (IPCC, 2014) influence the stability and biodiversity of ecological networks. Increasing the stability of food webs is possible via obtaining weaker interaction strengths (May, 1972) which is important for maintaining biodiversity. We found that loss of habitat complexity would lead to increased interaction strength via decreasing half-saturation density. Additionally, even though the constant interaction strengths through different patch sizes lead to constant population densities, it would result in a less absolute number of individuals in smaller patches. When there are a reduced number of individuals in smaller patches, decreasing patch size would cause species extinctions, especially at higher trophic levels. Therefore, shrinking patch sizes and homogenizing habitats would both lead to destabilization of ecological networks and biodiversity loss. Altogether, our study underlines the urgent need for protecting large complex habitats to save biodiversity.

Supplemental Information

CPP code for the individual based model

Includes a header file, a functions file and a main file to reproduce the model described in the main text.

DOI: 10.7717/peerj.2993/supp-1


Contains the ODD and further information about parameters.

DOI: 10.7717/peerj.2993/supp-2