Habitat selection patterns are density dependent under the ideal free distribution.

Despite being widely used, habitat selection models are rarely reliable and informative when applied across different ecosystems or over time. One possible explanation is that habitat selection is context-dependent due to variation in consumer density and/or resource availability. The goal of this paper is to provide a general theoretical perspective on the contributory mechanisms of consumer and resource density-dependent habitat selection, as well as on our capacity to account for their effects. Towards this goal we revisit the ideal free distribution (IFD), where consumers are assumed to be omniscient, equally competitive and freely moving, and are hence expected to instantaneously distribute themselves across a heterogeneous landscape such that fitness is equalised across the population. Although these assumptions are clearly unrealistic to some degree, the simplicity of the structure in IFD provides a useful theoretical vantage point to help clarify our understanding of more complex spatial processes. Of equal importance, IFD assumptions are compatible with the assumptions underlying common habitat selection models. Here we show how a fitness-maximising space use model, based on IFD, gives rise to resource and consumer density-dependent shifts in consumer distribution, providing a mechanistic explanation for the context-dependent outcomes often reported in habitat selection analysis. Our model suggests that adaptive shifts in consumer distribution patterns would be expected to lead to nonlinear and often non-monotonic patterns of habitat selection. These results indicate that even under the simplest of assumptions about adaptive organismal behaviour, habitat selection strength should critically depend on system-wide characteristics. Clarifying the impact of adaptive behavioural responses may be pivotal in making meaningful ecological inferences about observed patterns of habitat selection and allow reliable transferability of habitat selection predictions across time and space.

INTRODUCTION

Habitat selection is central to the ecological mission of understanding the contributory causes of variation in organismal abundance and distribution across space and time (Buckley et al., 2010; Matthiopoulos et al., 2015; Morris, 2003; Osorio‐Olvera et al., 2019; Pulliam, 2000; Rosenzweig, 1981, 1991). In heterogeneous environments, the availability and distribution of different habitats, as well as the capacity of organisms to adjust their habitat selection patterns accordingly, can drive vital rates and population dynamics. This notion has received much theoretical attention (Boyce et al., 2016; Buckley et al., 2010; Fryxell et al., 2005; Křivan, 1997; Matthiopoulos et al., 2015; Morris, 2003; Osorio‐Olvera et al., 2019; Pulliam, 2000) and some empirical support (Betini et al., 2015; Boyce & McDonald, 1999; Fryxell et al., 2004; Matthiopoulos et al., 2019; Street et al., 2017; Weber et al., 2017). Statistical models based on patterns of habitat selection are moreover extensively used as predictive tools in applied ecology, particularly in forecasting population geographical distributions under various climate‐ and landscape‐change scenarios. However, it is becoming increasingly clear that observed habitat selection patterns vary considerably across individuals, geographical domains and years depending on both spatial and temporal grain and extent (Johnson, 1980; Kleiven et al., 2018; Prokopenko, Boyce, et al., 2017; Randin et al., 2006; Soininen et al., 2018). Such variation is poorly understood, despite its potential significance, compromising the utility of habitat selection models as both inferential and predictive tools.

Consumer habitat selection patterns often vary with consumer abundance (Morris, 1987; Rosenzweig, 1981, 1991). As consumer density increases, resources within preferred habitats become limiting due to exploitation by conspecifics, so it can be beneficial for consumers to become less selective in their habitat use (i.e. expend their space use to include habitats of lower quality; see further discussion in the next section). This relationship is termed ‘negative density‐dependent habitat selection’, and has been empirically documented in a variety of species and systems (Bledsoe & Ernest, 2019; Falcy, 2015; McLoughlin et al., 2010; Morris, 1988; O'Neil et al., 2020; Paterson & Blouin‐Demers, 2018; Rosenzweig, 1981, 1991; Shepherd & Litvak, 2004; Sundell et al., 2012; van Beest, McLoughlin, et al., 2014; van Beest, Uzal, et al., 2014). Nevertheless, consideration of population‐density effects is mostly absent in the empirical habitat selection literature.

Habitat selection patterns also vary with changes in the availability of suitable habitats and/or resource abundance within these habitats, leading to a so‐called ‘functional response in habitat selection’ (also known as ‘availability’ or ‘resource’ dependence; Beyer et al., 2013; Godvik et al., 2003; Holbrook et al., 2019; Losier et al., 2015; Matthiopoulos et al., 2011; Mauritzen et al., 2014; Mysterud & Ims, 1998; Prokopenko, Boyce, et al., 2017). Such reaction norms are thought to stem from a nonlinear relationship between landscape‐scale resource availability and its relative utility to the consumer (due to saturation, threshold or trade‐off effects), and are often assumed to be monotonic, an assumption that lacks formal theoretical grounding (Holbrook et al., 2019; Matthiopoulos et al., 2011). In fact, the mechanisms invoked in explaining availability dependencies are often vague with regard to the nature of the shift in ‘resource availability’, which may refer to either a shift in the relative frequency of habitat classes, or a shift in the intrinsic value of these habitat classes, with their relative frequency remaining unchanged (Duparc et al., 2019; Paolini et al., 2019; Gaudry et al., 2018).

It is thus widely acknowledged that habitat selection patterns vary with densities of both consumers and resources, yet a lack of explicit theoretical underpinning hinders our capacity to predict the effects of density dependencies (McLoughlin et al., 2010; van Beest et al., 2016; but see Matthiopoulos et al., 2011). Because of the potential impact on the observed spatial relationship between consumers and their resources, our inability to account for these density dependencies clouds our interpretation of empirically parameterised models of animal space‐use patterns, particularly the widely‐used Habitat Selection Functions (HSFs; Manly et al., 2002). Consequently, ecological insights gained from HSFs may be site‐ and time specific, hindering their utility for predicting outcomes in other ecosystems or even the same ecosystem at a later time (Beyer et al., 2010; Bledsoe & Ernest, 2019; Matthiopoulos et al., 2011; Paolini et al., 2019; Street et al., 2017; Zurell et al., 2018).

Our aim here is to show how habitat selection patterns might be affected by consumer and resource densities based on a simple model of optimal habitat use. We start by briefly discussing the classical ideal free distribution (IFD, Fretwell & Lucas, 1969) perspective on density‐dependent habitat selection. We then proceed to describe the formal link between habitat‐dependent animal density and habitat selection inference, based on an inhomogeneous Poisson point process. We believe this is an essential step in the translation of a theoretical perspective into potentially useful empirical predictions. Our models demonstrate that, under a wide range of conditions, optimal habitat selection behaviour would be expected to lead to consumer and resource density dependence in habitat selection patterns.

The ideal free distribution and density‐dependent habitat selection

Several earlier researchers (Matthiopoulos et al., 2015; Morris, 2003; Pimm & Rosenzweig, 1981; Van Der Meer & Ens, 1997), have pointed out that a useful theoretical perspective on spatial ecology can be gained using the well‐known concept of the IFD (Fretwell & Lucas, 1969). IFD consumers are omniscient, equally competitive and freely moving, and are hence expected to distribute themselves across a heterogeneous resource landscape such that fitness is equalised across the population (Křivan et al., 2008). While these assumptions are clearly unrealistic, the IFD is an evolutionary stable strategy (ESS; Křivan et al., 2008), and often identifies many of the same spatial patterns that could be expected to emerge in general from more complex, adaptive patterns of behavioural decision‐making (Fryxell & Lundberg, 1998; Matthiopoulos et al., 2015; Morris, 1988, 2003; Rosenzweig, 1981, 1991; Sutherland, 1983; Van Der Meer & Ens, 1997). Moreover, the central assumptions underlying IFD theory are fully compatible with the (often implicit) ‘pseudo‐equilibrium’ and ‘unrestricted‐availability’ assumptions underlying commonly used statistical models of species distribution and habitat selection (Guisan & Thuiller, 2005; Johnson, 1980; Lele et al., 2013; Matthiopoulos et al., 2015).

What might the IFD tell us about density‐dependent habitat selection? The classical IFD construct relies on viewing consumers as selecting among a set of available habitats that differ only in their intrinsic quality. As long as total consumer abundance is low, only high‐quality (resource‐rich) habitats should be occupied (‘used’), but as consumer abundance increases, the effects of competition on the rate of resource acquisition within occupied habitats, compromising fitness returns, would be expected to lead to expansion of the range of acceptable habitats. The IFD should accordingly result in positive density‐dependent habitat occupancy (a greater range of habitats are occupied at higher population density), and hence negative density‐dependent habitat selection (Fretwell & Lucas, 1969; Morris, 2003; Rosenzweig, 1981, 1991; Shepherd & Litvak, 2004). On the other hand, it is commonly assumed that relative differences in intrinsic quality among different habitats never vary, that is, the within‐habitat rate of density‐dependent decline in fitness is constant. Consequently, the relative density (distribution) of IFD consumers distributed among occupied habitats is expected to be independent of total consumer abundance (Falcy, 2015; Kacelnik et al., 1992; Lessells, 1995; McLoughlin et al., 2010; Morris, 1988; O'Neil et al., 2020). Hence, whether the classical IFD model results in density‐dependent habitat selection patterns rests on one's frame of reference, as well as on the definitions of ‘distribution’ and ‘selection’ (Matthiopoulos et al., 2015). Lastly, if the within‐habitat rate of density‐dependent decline in fitness is not constant, the relative density of IFD consumers within occupied habitats is expected to be density dependent (Cressman & Křivan, 2006; Křivan, 2003; Morris, 2003; Van Der Meer & Ens, 1997). In the following sections we will demonstrate that it is only under a very restrictive set of assumptions that we might expect the IFD to be independent of consumer and resource densities.

MATERIALS AND METHODS

Inhomogeneous Poisson SPP and habitat selection

Spatial point processes (SPP) are models describing the probability of a point event (e.g. animal presence) across space (Cressie, 2015). A spatial point pattern, such as the collection of points where a species occurs, is a stochastic realisation of a SPP. Under a homogeneous Poisson SPP, the probability of observing exactly N events within a spatial unit of size m is: , where λ is the (homogeneous) Poisson SPP intensity (the expected density of events anywhere within the spatial domain; Cressie, 2015). The inhomogeneous Poisson SPP is the simplest alternative to the complete spatial randomness assumed by the homogeneous Poisson SPP and is hence a fundamental model in spatial ecology (Cressie, 2015; Hooten et al., 2017). Under an inhomogeneous Poisson SPP, intensity varies across spatial domain Ω, so that the local SPP intensity in location x (ϵΩ) is:Here μ is a normalisation constant ensuring that the mean (homogeneous) intensity across Ω is . The β's are parameters linking the relative SPP intensity at location x to the local values of q habitat dimensions, H =1: (x), such as resource availability (e.g. nesting sites) or autecological conditions (e.g. solar radiation). The inhomogeneous Poisson SPP thus provides a crucial statistical link between observed spatial occurrence patterns and the underlying drivers of spatial heterogeneity in population density (i.e. habitat quality); for habitat unit of size m centred at x, λ(x) = E[N(x)]/m (Cressie, 2015; Hooten et al., 2017; Renner et al., 2015; Velázquez et al., 2016; Warton & Shepherd, 2010).

Habitats are sectors of environmental space differing along one or more of its dimensions, such as resource availability (e.g. food, water, or shelter) or autecological conditions (e.g. temperature, pH or salinity; the constituents of the species' fundamental niche), and hence in their quality (Aarts et al., 2008; Morris, 2003). Habitat units are discrete, environmentally homogeneous, sectors of geographical space. Habitat selection is the usage of a given habitat more often than would be expected based on the availability of its units (Avgar et al., 2017; Lele et al., 2013; Matthiopoulos et al., 2020). The exponential Habitat Selection Function (eHSF) is a common statistical model of animal space use, yielding values that are proportional to the usage of a habitat unit centred at x: (Manly et al., 2002). In the context of eHSF, β is termed the selection coefficient for H, and it is readily estimable using standard logistic regression with presence‐background (also known as used‐available /‐control /‐pseudoabsence) data. Although other HSF formulations exist, eHSF is by far the most popular among field biologists due to its ease of application and the prevalence of presence‐only and remote‐sensing data (Avgar et al., 2017; McDonald et al., 2013).

In recent years, several authors have shown that the use of logistic regression with presence‐background data to fit an eHSF corresponds to estimating an inhomogeneous Poisson SPP (Aarts et al., 2012; Fithian & Hastie, 2013; Johnson et al., 2013; McDonald et al., 2013; McDonald, 2013; Renner et al., 2015; Warton & Aarts, 2013). This can be demonstrated by expressing p(x is used) as the Poisson probability that at least one event occurs within a habitat unit of size m (centred at x): , providing a simple formal link between ‘abundance’ and ‘occupancy’. As long as (because either μ or m are very small), is well approximated by . Based on Equation 1, we can now write:The selection coefficients in an eHSF are thus asymptotically unbiased estimators of the parameters of an inhomogeneous Poisson SPP (assuming habitat availability is correctly defined; Johnson, 1980; Matthiopoulos et al., 2020).

Assuming an inhomogeneous Poisson SPP, and for any two habitats, i and j (each a unique location in environmental space with coordinates and ), with corresponding total areas M and M (the sum of the areas of their respective habitat units), we can express the ratio of expected habitat‐specific densities, ρ and ρ, using Equation 1:Assuming our two habitat units are identical along all other habitat dimensions but k, and taking the natural logarithm of both sides, we obtain the following equation:Equation 4 is the mathematical definition of the eHSF coefficient, β, also known as the log relative‐selection strength (Avgar et al., 2017); it is the natural logarithm of the inter‐habitat (expected) density ratio divided by the inter‐habitat distance along the focal habitat dimension, k. Note that, for categorical habitat partitions (e.g. H ,  = ‘forest’ and H ,  = ‘meadow’), the inter‐habitat distance in environmental space is always 1 (; where is an indicator function, accordingly valued 0 or 1). In other words, β is a measure of the proportional population density change across the k'th dimension of environmental space.

Assuming IFD, the expected density ratio, and hence β, can be explicitly calculated as the fitness‐equalising solution across a two‐habitat system (see next section). Hence, under the prevalent assumptions of an adequately defined availability domain (Johnson, 1980; Matthiopoulos et al., 2020), and log‐linear relationship between λ and H =1: (Renner et al., 2015), β is a simple function of the IFD; if the IFD is consumer‐ and/or resource density dependent, so is exponential habitat selection strength.

A multidimensional model of the ideal free distribution

Our approach is based on the premise that ‘habitat quality’ should be more generally regarded as a function of multiple factors (multiple dimensions in environmental space), with various degrees of density‐dependent payoffs (Lampert et al., 2003; Matthiopoulos et al., 2015; Tyler & Gilliam, 1995). We refer to this extension as the multidimensional Ideal‐Free Distribution (hereafter, mIFD) as it can be applied to arbitrarily complex habitat‐dependent payoff functions, in contrast to the common perspective of the IFD as emerging from a single density‐dependent relationship.

Our model is framed around an ideal free consumer whose fitness gain ω within a specific habitat i depends on the density of resources (consumables) R, the density of conspecific consumers ρ and autecological conditions (non‐consumables) A. We express the fitness gain of consumers at i in the following manner:Here is the direct (additive) fitness gain from local autecological conditions, whereas is the contribution to fitness obtained from acquisition and assimilation of resources, which is a function of habitat‐specific resource density, consumer density and a vector of parameters, θ, that may or may not vary with autecological conditions. These parameters might include, for example, habitat‐specific variation in consumer search rate or net‐energy gain from a single consumption event. Hence, to keep our formulation as general as possible, we allow autecological conditions to affect ω via two different pathways, a direct (additive; as in Matthiopoulos et al., 2015) pathway, and an indirect (multiplicative; as in Pulliam, 2000) pathway operating through autecological effects on the parameters of . For example, temperature or habitat structure (e.g. vegetation cover or ruggedness) may directly influence consumer fitness (by affecting metabolic rates, thermoregulation costs or predation risk), but could also have multiplicative influences by affecting the parameters of , such as search rate and intraspecific competition intensity (e.g. due to temperature‐ or habitat‐dependent movement costs).

The spatial distribution of consumers conforms to an IFD provided that (a) their fitness‐gains are equalised across all occupied habitats and (b) no habitat with higher fitness exists. For simplicity, we assume that no habitat with higher fitness exists and we focus on the simple ‘snapshot’ IFD perspective, where R is the temporally static ‘standing stock’ of resources and conspecifics in i (Van Der Meer & Ens, 1997). Assuming that fitness within any single habitat is defined by Equation 5, it accordingly follows that for any two occupied habitats, i and j, the following fitness‐balancing relationship would be expected under IFD:As long as consumers occupy both habitats (i.e. both ρ and ρ are positive), we can represent the resulting mIFD as the log ratio of consumer densities that satisfies Equation 6: ln[ρ] – ln[ρ]. If the absolute value of ln[ρ] – ln[ρ] shifts with the magnitude of either the system‐wide density of consumers ((ρ + ρ)/2) or resources ((R + R)/2), the IFD, and hence β (see Equation 4), are consumer‐ or resource density dependent.

To investigate the ecological conditions leading to density‐dependent mIFD and eHSF we must first define an explicit functional form for (Van Der Meer & Ens, 1997). We set to one of the two commonly studied consumer‐dependent Type‐II formulations; either a Beddington–DeAngelis functional response (Beddington, 1975; DeAngelis et al., 1975), or an Arditi–Akcakaya functional response ( Arditi & Akçakaya, 1990). Both versions have been well‐documented under experimental conditions (Kratina et al., 2009; Novak et al., 2017; Prokopenko, Turgeon, et al., 2017). Under the Beddington–DeAngelis model, resource consumption rate is given by:Under the Arditi–Akcakaya model, resource consumption rate is given by:In both formulations, s[A] is the habitat‐specific fitness gain from a single consumption event (assimilation efficiency), a[A] is the habitat‐specific search rate, b[A] is the habitat‐specific time it takes a consumer to handle or process one resource item (b[A] = 0 corresponds to a Holling Type I functional response) and c[A] is the habitat‐specific effect of interference/competition/facilitation on a given consumer. Both formulations reduce to the classical consumer‐independent functional response (Holling Type II) if c[A] = 0, and both imply among‐consumer facilitation if c[A] < 0 (two situations that we will not consider here due to their irrelevance to IFD; Van Der Meer & Ens, 1997). Lastly, if c[A] = 1, the Arditi–Akcakaya model becomes the purely ratio‐dependent model of Arditi and Ginzburg (1989).

Substituting Equations 7 or 8 into Equation 6, and solving for ρ given ρ, R, R, θ[A] − θ[A], and g[A] − g[A], allows us to explore the expected relationship between mIFD (habitat selection strength) and the system‐wide densities of consumers and resources. For the sake of simplicity and gerality, we avoid specifying functional forms for θ[A] and g[A]. Instead, we solve the mIFD for pairs of habitats (i and j) that differ by a single unit of one parameter or variable (a[A], b[A], c[A], s[A], R or g[A]; dimensions in environmental space), but are identical in all other aspects. Our analysis thus considers 12 ecological scenarios, six types of between‐habitat differences under each of the two functional response formulations. Given the simplest imaginable ecological setting of a single type of resource, two habitats may minimally differ by their associated resource abundance, search rates, handling times, resource assimilation efficiencies, intensities of consumer interference or autecological fitness gains. Our aim here is to provide a theoretical perspective that is ‘as simple as possible, but not simpler’ by providing an exhaustive account of all possible configurations of a two‐habitat system under a generic fitness function. In Appendix S1 we provide the explicit terms for ρ 1 for each of 12 scenarios; the six types of between‐habitat differences (a[A], b[A], c[A], s[A], R or g[A]) under each of the two functional response formulations considered here. In Figures 1, 2, 3 we plot the mIFD‐based habitat selection strength (β = ln[ρ] − ln[ρ]; Equation 4 where H ,  − H ,  = 1) as a function of mean consumer and resource densities ((ρ + ρ)/2 and (R + R)/2 respectively). As ours is a qualitative investigation, parameter values, as well as variable magnitudes, are arbitrary. In the plots shown below, unless stated otherwise, all parameters were set to 1 except for assimilation efficiency (s[A]), which was set to 10.

FIGURE 1

mIFD‐based habitat selection strength as function of mean consumer density (x‐axis) and mean resource density (dotted line: (R + R)/2 = 1, dashed line: (R + R)/2 = 5, solid line: (R + R)/2 = 10). The different sub‐plots represent different ecological scenarios: an Arditi–Akcakaya model where the two habitats differ by instantaneous resource availability only (R + R)/2 = 1; a), handling time only (b[A] = 1.5, b[A] = 0.5; c), or autecological fitness gain only (g[A] = 1.5, g[A] = 0.5; e), or a Beddington–DeAngelis model where the two habitats differ by instantaneous resource availability only (R − R = 1; b), handling time only (b[A] = 1.5, b[A] = 0.5; d), or autecological fitness gain only (g[A] = 0.5, g[A] = 0.5; f)

FIGURE 2

mIFD‐based habitat selection strength as function of mean consumer density (x‐axis), where the two habitats differ only by resource assimilation efficiency (s[A] – s[A] = 1). The different sub‐plots represent different ecological scenarios: a type II (b[A] = b[A] = 1) Arditi–Akcakaya model across three different mean resource densities (a; dotted line: R = R = 1, dashed line: R = R = 5, solid line: R = R = 10), a type II (b[A] = b[A] = 1) Beddington–DeAngelis model across three different resource densities (b; dotted line: R = R = 1, dashed line: R = R = 5, solid line: R = R = 10), a type I (b[A] = b[A] = 0) Arditi–Akcakaya model across three different magnitudes of assimilation‐efficiencies (c; thin line: s[A] = 1.5, medium line: s[A] = 10.5, wide line: s[A] = 100.5), and a type I (b[A] = b[A] = 0) Beddington–DeAngelis model across three different magnitudes of assimilation‐efficiencies (d; thin line: s[A] = 1.5, medium line: s[A] = 10.5, wide line: s[A] = 100.5)

FIGURE 3

mIFD‐based habitat selection strength as function of mean consumer density (x‐axis). The different sub‐plots represent different ecological scenarios: an Arditi–Akcakaya model where the two habitats differ by search rate only (a[A] − a[A] = 1) across three different magnitudes of search‐rate (a; thin line: a[A] = 1.5, medium line: a[A] = 5.5, wide line: a[A] = 10.5), a Beddington–DeAngelis model where the two habitats differ by search rate only (a[A] − a[A] = 1) across three different magnitudes of search‐rate (b; thin line: a[A] = 1.5, medium line: a[A] = 5.5, wide line: a[A] = 10.5), an Arditi–Akcakaya model where the two habitats differ by consumer interference only (c[A] − c[A] = 1) across three different magnitudes of consumer interference (c; thin line: c[A] = 1.5, medium line: c[A] = 5.5, wide line: c[A] = 10.5), and a Beddington–DeAngelis model where the two habitats differ by consumer interference only (c[A] − c[A] = 1) across three different magnitudes of consumer interference (d; thin line: c[A] = 1.5, medium line: c[A] = 5.5, wide line: c[A] = 10.5)

RESULTS

In classical ‘snapshot’ IFD models (no resource dynamics; no autecological differences; purely ratio‐dependent functional response), the IFD, ln[ρ] – ln[ρ], depends on the ratio of resource densities, ln[R] – ln[R], but not on the absolute densities of either consumers or their resources (Křivan, 2003; Lessells, 1995; Van Der Meer & Ens, 1997). Hence, even under this simplest of models, the IFD shifts with absolute resource density if the latter is not proportional to the inter‐habitat density difference (i.e. if (R – R)/(R + R) is not a constant for all R's and R's). Resource‐selection strength, β (; Equation 4), always shifts with absolute resource density (Figure 1a).

More generally, out of the 12 different scenarios considered here (six types of between‐habitat differences, and two functional responses), nine scenarios result in consumer density dependence, whereas eight scenarios result in resource density dependence in habitat selection strength (Figures 1, 2, 3). Only two scenarios result in no density dependence, an Arditi–Akcakaya model where the two habitats differ by search rate only (a[A] − a[A] = 1; Figure 3a), and a Beddington–DeAngelis model where the two habitats differ by consumer interference only (c[A] − c[A] = 1; Figure 3d). Note, however, that density‐independent habitat selection strength is also expected under a Type I Arditi–Akcakaya model (i.e. in the absence of handling time; b[A] − b[A] = 1) where the two habitats differ by resource assimilation efficiency only (s[A] − s[A] = 1; Figure 2c). That said, even in these density‐independent scenarios, habitat selection strength varies with the magnitude of the autecological effect (a[A] + a[A], c[A] + c[A] and s[A] + s[A]).

Across all density‐dependent scenarios, the response is always nonlinear, with habitat selection strength either exponentially increasing or decreasing with consumer density, resource density or both (Figures 1, 2, 3). In the case of the scenarios assuming inter‐habitat difference in autecological fitness gain (g[A] − g[A] = 1), the response is also non monotonic, first decreasing and then increasing with both consumer and resource densities (Figure 1e,f). The shapes of these response curves are qualitatively insensitive to the absolute magnitude of parameters or variables, or to whether or not the inter‐habitat difference is linearly proportional to the mean value (i.e. whether or not ). The only exception is a type I (b[A] − b[A] = 0) Arditi–Akcakaya model in the case of an inter‐habitat difference in autecological fitness gain, which results in a monotonic increase in β with consumer density and a monotonic decrease in β with resource density.

DISCUSSION

We have shown that the distribution of adaptive consumers should shift with the overall intensity of consumer competition (density dependence) and resource density (availability dependence). Such distributional shifts are manifested as density‐dependent shifts in the estimated coefficients of habitat selection, at least as it is inferred based on the inhomogeneous Poisson point process. These shifts are nonlinear and often non‐monotonic, making them potentially challenging to interpret when dealing with empirical data. Furthermore, as empirical applications typically examine distributional patterns across vast landscapes, which may include habitats that are either unoccupied and/or function as sinks (i.e. have non‐positive fitness), density‐dependent distributional shifts may prove even more extreme than depicted here. To the extent that organisms in real systems exhibit adaptive spacing patterns across the landscape that they occupy, our findings suggest that spatial patterning in consumer‐resource systems should be conditional on both spatial habitat heterogeneity and changes over time or space in mean consumer or resource abundance. Consequently, it may be extremely challenging to interpret and predict ecological requirements, habitat suitability and the potential effects of habitat heterogeneity on population interactions based on simple snapshots of spatial distribution patterns.

Our findings reinforce earlier conclusions drawn from incorporating consumer interference into IFD (Van Der Meer & Ens, 1997) and from isodar theory (Morris, 2003). Isodars are lines of equal fitness across the state‐space of two alternative habitats (i.e. population density in one habitat plotted against population density in the other), a graphical manifestation of density‐dependent habitat selection (Morris, 1988, 1989, 1990). Despite its relevance to our understanding of density‐dependent, single‐species, habitat selection patterns (Morris, 2003), isodar theory was developed and employed primarily in the context of scale‐dependent habitat distribution (Morris, 1992), and multi‐species coexistence/competition (Guthrie & Moorhead, 2002; Tarjuelo et al., 2017). Nevertheless, in the context of single‐species habitat selection, ‘Isodars need not be linear. Interference, for example, can result in curved fitness functions. … (1) multiple strategies of habitat distribution will exist whenever isodars are curved, or when they are linear with non‐zero intercepts…’ (Morris, 2003). Such curved isodars would be expected whenever density dependence is nonlinear, such as in the functional response models used here, in the Gilpin and Ayala (1973) general‐logistic model used in Morris (2003), or in any of the interference models used in Van Der Meer and Ens (1997). Isodars with non‐zero intercepts are expected whenever habitats differ in their fitness returns at very low densities, such as when habitats differ in autecological (density‐independent) conditions. These parallels highlight the central message of both our current work and these previous theoretical contributions: under the typical ecological conditions of nonlinear functional responses and consumer interference, habitat selection is expected to be a non‐linear and often non‐monotonic function of mean consumer population density and/or mean resource density.

The most important consequence of this fundamental mathematical reality is that there is no reason to expect exponential habitat selection models to be transferable across systems differing in mean consumer density or mean resource density. At the very least, this suggests that habitat selection and species distribution studies might benefit from critical evaluation of model transferability, particularly when the research objective is ecological forecasting (Radchuk et al., 2019). Fortunately, critical evaluation techniques are readily available (Bahn & Mcgill, 2013; Fieberg et al., 2018; Journé et al., 2019; Kleiven et al., 2018; Randin et al., 2006; Soininen et al., 2018; Yates et al., 2018), and should be used more often.

The multidimensional IFD model we have outlined here is compatible with previous IFD‐based models and analysis, suggesting that selective space use is not only a function of the quality of a single habitat attribute, but of the integrated quality and spatial configuration of all available habitats (Abrahams & Dill, 1989; Morris, 2003). A simplified version of our multidimensional IFD model, termed ‘IFD with costs’, was proposed and experimentally evaluated more than 20 years ago by Tyler and Gilliam (1995). The concept was then further developed and evaluated by Lampert et al. (2003), who also showed that Daphnia population distribution is skewed less towards food‐rich areas in experimental mesocosms when the temperature is lowest in depth zones with highest food abundance (i.e. there is an inverse correlation between resource abundance and autecological conditions favouring fitness). Lampert (2005) then elegantly demonstrated how IFD with costs translated into density‐dependent spatial distribution in Daphnia. Similarly, Halliday and Blouin‐Demers (2014) experimentally demonstrated that flour beetles adaptively balance thermoregulation and resource acquisition via density‐dependent habitat selection. Finally, Matthiopoulos et al. (2015) presented a comprehensive theoretical framework establishing the link between IFD‐based habitat selection and animal population dynamics across a multidimensional environmental space.

These theoretical and experimental studies reinforce a growing body of observational studies suggesting that habitat selection patterns should be expected to vary with consumer and resource densities, and thus be context‐dependent, with obvious ramifications for what can and cannot be learned from habitat selection models. Consider the common ecological scenario where consumers are faced with a spatial trade‐off between safety and food resources. Because both per capita predation risk and per capita rates of resource consumption typically decrease with consumer abundance, consumers are expected to shift their selection from safety to habitats with high food abundance as their population increases in size. Inference based on observations of ecosystems with low consumer density may thus be highly misleading when applied to ecosystems with a high density of consumers, and vice‐versa. Furthermore, understanding the influence of variation in consumer and resource densities on habitat selection could play a key role in understanding consumer‐resource dynamics, as it defines how per capita consumption rate will vary over space and time (Matthiopoulos et al., 2015; Schmitz et al., 2017). Although we have limited our derivation to the simple ‘snapshot’ model where resources undergo neither depletion nor growth, we note that once resource dynamics are explicitly considered, the IFD may be either density dependent or density independent, depending on the details of these dynamics (Křivan, 2003). It could prove useful (albeit challenging) to extend our model into a coupled dynamical system where both consumers and their resources are distributed in an ideal‐free manner, or where a given species is both a resource (for a predator or pathogen) as well as a consumer, such as the case for nearly all herbivores (Abrams, 2007; Křivan, 1997, 2003; Křivan et al., 2008; Lessells, 1995).

While our assessment of the lack of transferability of habitat selection models may seem a little depressing, our model provides a theoretically grounded starting point that could be helpful in developing a set of truly robust predictions based on more realistic models for the kinds of context‐dependent habitat selection that are often observed in real ecosystems (Beyer et al., 2013; Bledsoe & Ernest, 2019; Falcy, 2015; Godvik et al., 2003; Holbrook et al., 2019; Losier et al., 2015; Matthiopoulos et al., 2011; McLoughlin et al., 2010; Paterson & Blouin‐Demers, 2018; Prokopenko, Boyce, et al., 2017; Shepherd & Litvak, 2004; Sundell et al., 2012; van Beest et al., 2016; van Beest, Uzal, et al., 2014). IFD models represent one of the simplest ways of thinking about adaptive space use, because they are based on several unrealistic assumptions. Simple models often have the advantage, however, of focusing our attention on key properties that are upheld in more complex and realistic model variants. First, consider our assumption that consumers are both ideal and free, clearly a debatable assumption at best. It is nonetheless difficult to see how relaxing this assumption to accommodate more realistic behaviour would counteract our qualitative predictions. For example, the added behavioural complexity of territoriality should result in even stronger density dependencies (O'Neil et al., 2020). Secondly, by considering only unidimensional differences between habitats, we have assumed that there is no spatial correlation across different habitat dimensions (e.g. that resource rich habitats are not consistently better or worse autecologically). In reality, many ecological landscapes will be characterised by multiple positive or negative cross‐correlations. If the spatial correlation between different fitness correlates is negative, fitness‐optimising consumers face a trade‐off that depends on the relative magnitude of their selection coefficients. Since we have shown here that these coefficients are context dependent, but not in the same fashion, consumers might appear to avoid otherwise favourable autecological conditions in some contexts, and abundant resources in others. In other words, in the presence of spatial correlation across habitats, the mechanism described here could very well lead to qualitative shifts in the inferred direction of habitat selection patterns.

A potential step forward in overcoming the challenges presented by theses context dependencies in habitat selection patterns, is to invest more in obtaining estimates of system‐wide population size (of the focal species as well as of other interacting species) and resource availability. Measuring variability in these state variables requires replicated distribution studies across space and/or time. Emerging remote‐sensing technologies, such as wildlife cameras and drones, could enhance our capacity to obtain such complex data. Once obtained, these estimates could be incorporated into more appropriate forms of statistical analysis where context dependency could be adequately controlled for (Holbrook et al., 2019; Matthiopoulos et al., 2011; O'Neil et al., 2020). On the other hand, if empirical patterns of density‐dependent habitat selection prove as extreme and variable as the theoretical patterns shown here, statistical adjustments alone may fall short of improving model transferability (Radchuk et al., 2019; Yates et al., 2018). It may thus be necessary to employ more mechanistic approaches, either explicitly modelling the link between population distribution and fitness (Buckley et al., 2010; Matthiopoulos et al., 2015, 2019; Osorio‐Olvera et al., 2019; Pulliam, 2000), or explicitly modelling movement behaviour by individuals in response to fitness variation across ecologically meaningful landscapes (Avgar et al., 2013, 2015, 2016; Fryxell et al., 2020; Jonsen et al., 2018; Patterson et al., 2008, 2017). Whether such mechanistic approaches could aid in obtaining robust inference on animal space‐use patterns and generating transferable predictive models is still very much an open question. Further research is needed to, first, establish the theoretical conditions under which it is at all possible to forecast space‐use patterns, and second, to empirically evaluate habitat selection model transferability in real ecological systems.

AUTHORS' CONTRIBUTIONS

T.A. designed the model and wrote the manuscript with much input from J.M.F. and G.S.B.

Appendix S1

Click here for additional data file.