Climate change and integrodifference equations in a stochastic environment
Climate change impacts population distributions, forcing some species to migrate poleward if they are to survive and keep up with the suitable habitat that is shifting with the temperature isoclines. Previous studies have analyzed whether populations have the capacity to keep up with shifting temperature isoclines, and have mathematically determined the combination of growth and dispersal that is needed to achieve this. However, the rate of isocline movement can be highly variable, with much uncertainty associated with yearly shifts. The same is true for population growth rates. Growth rates can be variable and uncertain, even within suitable habitats for growth. In this paper we reanalyze the question of population persistence in the context of the uncertainty and variability in isocline shifts and rates of growth. Specifically, we employ a stochastic integrodifference equation model on a patch of suitable habitat that shifts poleward at a random rate. We derive a metric describing the asymptotic growth rate of the linearized operator of the stochastic model. This metric yields a threshold criterion for population persistence. We demonstrate that the variability in the yearly shift and in the growth rate has a significant negative effect on the persistence in the sense that it decreases the threshold criterion for population persistence. Mathematically, we show how the persistence metric can be connected to the principal eigenvalue problem for a related integral operator, at least for the case where isocline shifting speed is deterministic. Analysis of dynamics for the case where the dispersal kernel is Gaussian leads to the existence of a critical shifting speed, above which the population will go extinct, and below which the population will persist. This leads to clear bounds on rate of environmental change if the population is to persist. Finally we illustrate our different results for butterfly population using numerical simulations and demonstrate how increased variances in isocline shifts and growth rates translate into decreased likelihoods of persistence.
Keywords:integrodifference equations climate change, forced migration speed stochastic environment population persistence
Msc:45R05 45C05 92D25
The consequences of climate change on population abundance and distribution have been widely investigated for the last two decades. One of these consequences is a modification of range distributions. Indeed, we know that, for a large variety of vertebrate and invertebrate species, climate change induces range shift toward the poles or higher altitudes, contraction or expansion of the habitat and habitat loss (Parmesan and Yohe 2003; Hickling et al 2006; Menendez et al 2014; Parmesan 2006; Lenoir et al 2008, amongst other). Mathematical models and simulations that include climate change have predicted an effect of climate change on range distribution through habitat migration, habitat reduction and expansion and habitat loss (Parmesan and Yohe, 2003; Ni, 2000; Hu et al, 2015; Malcolm and Markham, 2000; Polovina et al, 2011; Parr et al, 2012; Hazen et al, 2013). In this paper we are interested in understanding, with the aid of a mechanistic model, the effect of shifting range on population persistence when the yearly range shifts and the population growth rates are stochastic.
Mechanistic models have been used to study the effect of shifting boundaries on population persistence, by considering a suitable habitat, which is a bounded domain where the population can grow, that is shifted toward the pole at a forced speed . Potapov and Lewis (2004) used a reaction-diffusion system with a moving suitable habitat to investigate the effect of climate change and shifting boundaries on population persistence when two populations compete with one another. Berestycki et al (2009) used a similar equation, in a scalar framework to study the persistence property of one population facing shifting range, and characterised persistence as it depends on the shifting speed. More recent papers also investigate the effect of shifting range on population persistence of single populations (Richter et al, 2012; Leroux et al, 2013; Li et al, 2014). In all these reaction-diffusion equations are used to model the temporal evolution of the density of a population in space and time. That is, individuals are assumed to disperse and grow simultaneously. In these models dispersal is local in the sense that the population disperses to its closest neighbourhood, in a diffusive manner.
Another approach to modelling the temporal evolution of the density of a population, is to consider populations that disperse and reproduce successively, and to allow for nonlocal dispersal. In this case, integrodifference equations are the appropriate model for the dynamics of the density . integrodifference equations, introduced by Kot and Schaffer (1986) to model discrete-time growth-dispersal, assume time () to be discrete, and space to be continuous. From one generation to the next the population grows, according to a nonlinear growth function and then disperses according to a dispersal kernel , so
Strictly, the dispersal kernel is a probability density function describing the chance of dispersal from to .
We consider a self-regulating population, with negative density-dependence so the slope of the growth is assumed to be monotically decreasing with for and for , where is the carrying capacity. As we consider a population that is not subject to an Allee effect, the standard assumption on the growth function is that the geometric growth rate is the largest at lowest density, that is achieves its supremum as approaches 0. We denote . When we wish to explicitly distinguish between populations with different geometric growth rates we modify our notation, replacing by . Within this framework, we consider two types of growth dynamics: compensatory and over-compensatory. The compensatory growth dynamics are monotonic with respect to density whereas the overcompensatory growth dynamics have a characteristic “hump” shape. (Figure 1).
The eigenvalue problem associated with the linearisation of (1.1) about is
and persistence of the population depends upon whether falls above or below one (Kot and Schaffer, 1986; VanKirk and Lewis, 1997; Lutscher and Lewis, 2004). An approximate method for calculating employs the so-called dispersal success approximation (VanKirk and Lewis, 1997). Without loss of generality one can assume that and integrating the previous equation on we get
where is the dispersal success. This function represents the probability for an individual located at to disperse to a point within the domain. The so-called dispersal success approximation thereby allowing the principal eigenvalue to be be estimated by
(VanKirk and Lewis, 1997). This approximation gives as the growth rated times the estimated proportion of individuals that stay within the suitable habitat from one generation to the next. A modified dispersal success approximation recently introduced by Reimer et al (2016) improves upon the dispersal success approximation assumption that the population is uniformly distributed within the favourable environment. Reimer et al (2016) introduced a modified approximation that weights the dispersal success values by the proportion of the population at each point. They defined the modified dispersal success approximation by
and showed that this gave a better approximation to the eigenvalue. We will employ both versions of the dispersal success approximation (1.4) and the modified dispersal success approximation (1.5) in our calculations later in this paper.
So far we have considered only the dependence of the growth on the density . In the framework for climate change, the population can grow differently depending on where it is located with respect to space and time. To take this into account we introduce a suitability function, (), which depends on space and time and multiplies the growth map . As we consider populations whose suitable habitat shifts toward the pole, we choose a particular form for suitability function, , where is the initial suitability function in the absence of climate change and is a parameter standing for the center of the suitable habitat. The simple case where the habitat shifts at a constant speed is given by . However, as we describe below, it is also possible to allow to vary randomly about .
Following the growth stage, the population disperses in space according to the dispersal kernel . The kernel is assumed to be positive everywhere: that is, the probability of dispersing to any point in space is always positive. Dispersal kernels are typically assumed to depend only on the signed distance between two points in space, i.e. only on the dispersal location relative to the source location. If the population has no preferred direction of dispersal, the kernel is symmetric, depending only upon the distance between source and dispersal locations. This is not the case, however, in rivers where the population is subjected to a stream flow for example, or in environments with a prevailing wind direction that affects dispersal. In this paper we consider the Gaussian and the Laplace as examples of typical symmetric dispersal kernels (Figure 2)
In this framework of discrete-time growth-dispersal models, Zhou and Kot (2011) investigated the effect of climate change and shifting range on the persistence of the population and highlighted the possible existence of a critical shifting speed for persistence. More recently several works investigated the effect of shifting range in more general discrete time growth-dispersal models (Zhou and Kot, 2013; Harsch et al, 2014; Phillips and Kot, 2015).
Recent research reports an increase in the environmental stochasticity in population dynamics, partly due to climate change and its effect on the frequency and the intensity of extreme climatic events covering large areas of the globe (Saltz et al, 2006; IPCC, 2007; Kreyling et al, 2011). It is also known that the projected consequences for population ranges vary, depending on the different scenarios related to climate change (see, for example, IPCC (2014)).
In this paper we focus on the effect of environmental stochasticity on population persistence in the presence of a shifting range, using a stochastic population model. The development of stochastic population models in population ecology was initially motivated by the study of the effect of environmental stochasticity on population dynamics and on the large time behaviour of the population (May, 1973; Turelli, 1977).
We incorporate stochasticity into our modelling framework in two ways, first with respect to the shifting speed for suitable habitat and second with respect to the growth rate within the suitable habitat. The model for the shifting speed gives the center of the suitable habitat as , where is a random variable. Here we assume that is unknown but fixed, depending on the scenario considered for the severity of global warming () (Figure 3). Variability of the growth dynamics at each generation, for example due to weather conditions or extreme climate events, is included through stochasticity in the growth function . Our approach is to incorporate the randomness into the geometric growth rate (Figure 4) and so in any given year the growth function is given by where .
Our goal is to mathematically analyze the effect of climate change and shifting range on population dynamics for a species that grows and then disperses at each generation, taking into account stochasticity induced by environmental variability and climate change. In Section 2 we derive the model and state model assumptions. In Section 3 we first define a persistence condition, derived from the papers of Hardin et al (1988) and Jacobsen et al (2015). We then detail computation of the persistence criterion and highlight the link between persistence of the population and principal eigenvalue of the operator linearised about 0 for the case where the shifting speed is not random (). We also prove that, when the dispersal kernel is Gaussian and the speed is not random, there exists a critical shifting speed characterising persistence of the population. Speeds above this critical value will drive the population to extinction, while speeds below will allow the population to persist. In section 4 we apply the theory to an example in butterfly population subject to changing temperatures in Canada (Leroux et al, 2013) and use numerical simulations to investigate the dependence of the critical shifting speed on the variance of the dispersal Gaussian kernel. We also compute the persistence criterion as a function of the variance of the dispersal kernel for Gaussian and Laplace kernel. Lastly we numerically investigate the effect of the variability in the yearly shift and in the growth rate on the persistence of the population. In the Appendix we draw on classical theory of spreading speeds for stochastic integrodifference equations to aid development of a heuristic link between the speed of the stochastic wave in the homogeneous framework and the critical domain size and persistence condition.
2 The model
In this section we derive the mechanistic model used to study population persistence facing global warming and habitat shifts. We explain the different assumptions made for each component throughout the derivation of the model and conclude by explaining how it results in the problem in a moving environment.
As already stated in the introduction, we use the theory of integrodifference equations to model the temporal dynamics of the density of the population , as introduced by Kot and Schaffer (1986). In the classical homogeneous case, the model is as given in equation (1.1) with given by and with given, nonnegative, compactly supported and bounded.
We make the following assumptions regarding the dispersal kernel
For all , , .
is well defined, continuous, uniformly bounded and positive in .
The first hypothesis means that takes the form of a difference kernel and depends only on the signed distance between and . The second holds for typical kernel, such as the Laplace and Gaussian (Figure 2).
To include the effect of climate change on range distribution in this model, we assume that the suitability of the environment is heterogeneous in the sense that
where the function stands for the suitability of the environment at generation and location .
We make the following assumptions for the suitability function :
Denoting as the reference point on a suitable habitat, we assume that
is compactly supported, nonnegative, bounded by 1 and is non trivial in .
Biologically these assumptions mean that the suitable environment has a constant profile that is shifted by at generation .
The model then becomes
We denote as the support of , i.e . Notice that only the population located in contributes to the growth from generation to generation . To introduce environmental stochasticity in our model, we assume that is a random process and with a sequence of random functions. We interpret as the geometric growth rate of the population at low density. These two forms of environmental stochasticity emphasise the dependence of the range shift and the growth rate on the strength of the climate change. Moreover we can be more precise about the form of the shift variable . Indeed, we assume that for all , where is a constant representing the asymptotic shifting speed and is a random variable representing the environmental stochasticity in the shift from one year to the next. In the introduction we stated that the asymptotic shifting speed itself may be uncertain. However from now on we consider it to be a fixed constant and study the problem of persistence of the population for different possible values of the asymptotic shifting speed .
Denoting by and the set of possible outcomes for at each generation, we assume that the elements are independent, identically distributed and bounded by appropriate values, namely:
is a sequence of independent, identically distributed random variables, with distribution
There exists such that for all with probability 1
There exists such that for all with probability 1
We consider a self-regulating population, with negative density-dependence and make the following assumptions on the growth function :
For any such that
is continuous, with for all ,
There exists a constant such that, for all ,
for all positive continuous function ,
If constants such that then ,
is right differentiable at 0, uniformly with respect to .
We denote as the right derivative of at 0 and assume for now that
, from Hypothesis 3(iii)
, with .
Hypothesis 4(i) means that the population does not grow when no individuals are present in the environment. We also assume that the growth of the population is bounded (Hypothesis 4(ii)b.) and consider a population not subject to an Allee effect, that is, the geometric growth rate is decreasing with the density of the population (Hypothesis 4(ii)c.). Hypothesis 4(ii)d. holds for typical growth functions as the ones illustrated by Figure 1. Notice that because of Hypothesis 4(ii)a., we do not yet consider models with over-compensatory competition here. Indeed we will need this monotonicity assumption to prove large time convergence of the solution of our problem. Nevertheless, we can extend our result to the case with over-compensation where is not assumed to be nondecreasing and so Hypothesis 4(ii)a. no longer holds. More details are given in Corollary 1, stated in the next section. Hypothesis 4(iii)a. assumes that the geometric growth rate at zero density is bounded from above and below by some positive constants, while Hypothesis 4(iii)b. implies that the growth term at the maximal density is positive, for all the possible environments.
Finally, denoting by , we make a last assumption:
There exists such that
for all , nonnegative continuous function.
This assumption means that there exists an environment that has the better outcome than all other environments in terms of population growth.
The general problem becomes
We are interested in the large time behaviour of the density . Using a similar approach to the one of Zhou and Kot (2011), we would like to study the problem in the shifted environment, so as to track the population. Considering the moving variables , and letting be the associated density in the moving frame (where the reference point is the centre of the suitable environment that is shifted by at time ), we obtain a problem where the space variables and are also random variables. To simplify the analysis, we consider the problem in the “asymptotic” moving frame. That is we consider and and using Hypothesis 1(i), equation (2.4) can be written as
Letting be the associated density in the moving frame, we have
Dropping the bar, we obtain the following problem, in the moving frame,
where is given as a nonnegative, non trivial and bounded function. Moreover, as stated in Hypothesis 3(ii), for all , and thus defining
we only have to study equation (2.7) for all , for all . The problem in the moving frame becomes
where is given as a nonnegative, non trivial and bounded function. Notice that this problem is now defined on a compact set .
Many results already exist for the deterministic version of (2.9). When and is deterministic, depending only on , it has been shown that, if , the magnitude of the largest eigenvalue of the linearised operator around 0 determines the stability of the trivial solution and thus the persistence of the population (see Zhou and Kot (2011, 2013) for the one dimensional problem, Phillips and Kot (2015) for the two dimensional problem). Harsch et al (2014) study a similar deterministic integrodifference equation modified to include age and stage structured population and show that if the magnitude of the principal eigenvalue of the linearised problem around 0 exceeds one then 0 becomes stable and the population goes extinct. Zhou also studied the associated deterministic problem (2.7), for more general deterministic functions and (Zhou, 2013, chapter 4).
In this paper we are interested in similar questions about persistence of the population, but now for a shifting environment that moves at a random speed and a population that reproduces at a growth rate, chosen randomly at each generation.
Mathematically the assumptions on , and that we make are similar than those made by Jacobsen et al (2015) in their study of the question of persistence of a population in temporally varying river environment. They assumed that the dispersal kernel and growth terms are randomly distributed at each time step but also assumed asymmetric dispersal kernels to take the effects of water flow into account. They used the theory of Hardin et al (1988) to develop a persistence criterion for the general model
where , and are independently identically distributed random variables and related this persistence criterion to the long-term growth rate. In what follows, we use similar theory to study our integrodifference equation (2.9) where .
3 Persistence condition for random environments with shifted kernel
In this section we derive a criterion that separates persistence from extinction for a population facing random environments and climate change. We highlight the dependence of this criterion on the asymptotic geometric growth rate at low density on one hand and on the shifted dispersal success function on the other hand. We then develop a connection between persistence and the magnitude of the principal eigenvalue of a linearised operator for the case when only the growth rate is stochastic. This allows us to conclude on the existence of a critical shifting speed for the case of Gaussian dispersal kernel.
3.1 Conditions for persistence
We first describe the derivation of the deterministic criterion characterising persistence from extinction of the population.
Define to be the solution of problem (2.9). Then , the sequence of population density at each generation , is a random process and for each positive, takes values in , the set of continuous, nonnegative function defined on . Using the same notation as Hardin et al (1988), we can rewrite equation (2.9) as follows
for all .
We have the following theorem about the large time behaviour of the solution
Assume that , , and satisfy all the assumptions stated in Section 2 (Hypotheses 1 - 5), let be the solution of problem (2.9), with a bounded initial condition . Then converges in distribution to a random variable as time goes to infinity, independently of the initial condition , and is a stationary solution of (2.9), in the sense that
with a random variable taking its values in with distribution (Hypothesis 3(i)). Denoting by the stationary distribution associated with and the probability that , we have that or
where is the solution of the linearised problem around 0, i.e for all
The metric can be interpreted as the growth rate of the linearised operator up to time and its limit, representing the asymptotic growth rate of the linearised operator. We have the following theorem
Let be defined as in (3.3), then
If , the population will go extinct, in the sense that ,
If , the population will persist, in the sense that .
The proof of this theorem follows from Hardin et al (1988, Theorem 5.3) and Jacobsen et al (2015, Theorem 2), and it is explained in the Appendix A. From these two theorems one can then deduce a corollary, for the case where the growth function is not assumed to be monotonic anymore.
Let , , and satisfy the assumptions from Section 2 except assumption 4(ii)a, that is, is not assumed to be nondecreasing. Let be the solution of (2.9), with a bounded initial condition . Let be defined as in (3.3), then
If , the population will go extinct, in the sense that for all with probability 1,
If , the population will persist, in the sense that with probability 1.
First note that
where is defined in Hypothesis 4(ii)b. This implies that
for all , for all . Then define lower and upper nondecreasing functions that satisfies Hypotheses 4
so their slopes match that of at zero
and they satisfy the inequalities
for all . For instance, we can choose the nondecreasing function
Then denoting , respectively , as the solution of problem (2.9), with instead of , respectively instead of , such that , we have
If , then the population associated with the process goes extinct in the sense that , which implies, using again (3.11), that , for all , with probability 1.
This completes the proof of Corollary 1.
This means that the results about persistence and extinction extend to the case of over-compensatory growth as illustrated in Figure 1, even though we cannot draw conclusions about the large time convergence of the solution.
Note that the persistence criterion as given by Theorem 3.2 or Corollary 1 is difficult to analyse further. It requires the calculation of the asymptotic growth rate of the linearised operator under random environmental conditions (equation (3.3)). However we can proceed heuristically regarding a necessary condition for persistence of the population by comparing the asymptotic shifting speed with the asymptotic spreading speed of the population in a homogeneous environment. Consider the population dynamics in a homogeneous environment (). Equation (3.4) becomes
If we consider exponentially decaying solution of the form , , substitution into (3.12) yields
where is the moment generating function of
The pointwise growth rate of the solution is given asymptotically by
and is positive if where
We interpret as the spreading speed of a population undergoing dispersal and stochastic growth in a spatially homogeneous environment (Appendix B.1). In the case where is diminished, so it is not equal to one at all point in space, may no longer suffice to give a positive pointwise growth rate. On the other hand, if we consider a bounded, compactly supported initial condition, , and thus for all there exists such that , then the population , solution of (3.12) starting with the initial condition ,
The pointwise growth rate of the solution will thus be smaller that the pointwise growth rate of the solution and is necessary to give a positive pointwise growth rate of .
3.2 Computation of the persistence criterion
In this section, we detail possible computation of the persistence criterion and highlight the link between and the principal eigenvalue of the linearised operator around zero.
Using the analysis of the previous section one can determine whether a population can track its favorable environment and thus persist ( in Theorem 3.2 and Corollary 1) or cannot keep pace with the shifting environment and goes extinct ( in Theorem 3.2 and Corollary 1). It would be interesting to understand the dynamics of in terms of the different parameters of the problem. Using the definition of (3.3), we have
where , are the realisations for the geometric growth rate and the yearly shift at generation 0. Continuing,
where , are the realisations for the geometric growth rate and the yearly shift at generation 1. In a similar manner, one can then derive the more general formula
the geometric mean of the geometric growth rate at zero and
and using the strong law of large number we have
with probability 1 and thus obtain
This formula for highlights the dependence of the persistence criterion on the different parameters of the problem. Notice that the distribution of the growth rate at zero affects the first part whereas the distribution of the yearly shifts affects the second part of the formula. We further analyse the effect of the variation in or the variation in using numerical simulation in Section 4.
Persistence criterion for deterministic shifting speed
We turn our attention to further analysis of the persistence criterion , as given in (3.29), when we assume that the shifting speed is not random. That is, we assume that and thus randomness comes only from the growth term. We still consider a population that sees its favorable environment shifted at a speed , but this speed is not assumed to have random variation from one generation to the next. In this case and
Define the linear operator such that
As is positive and is compact, the operator is compact (Krasnosel’skii, 1964) and strongly positive, that is for any function , there exists such that for all . Then applying the Krein-Rutman theorem, it follows that this operator possesses a principal eigenvalue such that for all other eigenvalues and is the only eigenvalue associated with a positive principal eigenfunction . That is and satisfy
for all . One can always normalise such that . Now choosing the initial condition for equation (2.9) so that , we have
with defined in (3.28).
One has to approximate the principal eigenvalue to be able to conclude about the persistence of the population. One can refer to the paper by Kot and Phillips (2015) for the description and implementation of different numerical or analytical methods to compute the principal eigenvalue of a linear operator.
In the specific case of Gaussian kernel with variance , that is
the principal eigenvalue of the shifted linear operator can be expressed as a decreasing function of . Indeed, we know that and satisfy for all
The last equality is equivalent to
Thus one has that
where is the principal eigenvalue of the linear operator (defined by (3.31) when ) and
is decreasing with . Assuming that the parameters of the problem are chosen so that when the population persists, that is
We can thus state the following proposition
Let , and satisfy the assumption from Section 2. Assume also that and is a Gaussian kernel with variance , that is
There exists such that
The critical shifting speed for persistence thus increases with the expected geometric growth rate. Heuristically one would imagine that when is small enough relative to then should increase with the variance of the dispersal kernel, as the population becomes more mobile. On the other hand when becomes too large relative to then the population would disperse outside its favorable environment and then should decrease with the variance of the dispersal kernel.
one can also approximate the principal eigenvalue using the dispersal success approximation
as defined in (1.4) or the modified dispersal success approximation
as defined in (1.5), with the dispersal success function when . We obtained similar conclusion than Reimer et al (2016) when comparing the principal eigenvalue , its dispersal success approximation and its modified dispersal success approximation (Figure 5), observing that the modified dispersal approximation gives a more accurate estimate of the principal eigenvalue than the regular dispersal success approximation.
In this section we derived a analytical tool to characterise the persistence of a population facing shifting range in a stochastic environment. We also analysed further the dependence of the criterion with respect to the parameters of the model and investigated the existence of a critical shifting speed , separating persistence from extinction, in the particular case of Gaussian dispersal kernel.
4 Numerical calculation of persistence conditions, with application to shifting butterfly populations
We now focus on applying the stochastic model and theory to butterfly populations responding to climatic change. Our goal is twofold, first to illustrate the calculation of (3.29) and the critical shifting speed (3.46) in a ecologically realistic problem and second to investigate the impact of variability in environmental shift () and growth rate () on population persistence.
We assume that for each generation, there are only two possible environments, a good environment when and a bad environment when , such that
We consider , and with
To define good and bad environments we assume that the larger , the better the environment and the smaller the better the environment and thus consider
Leroux et al (2013) study the effect of climate change and range migration for twelve species of butterflies in Canada. In their paper the authors use a reaction-diffusion model to study the ecological dynamics of the population and its persistence properties. In their analysis they estimate the shifting speed and the growth rate for each species to compute the critical diffusion coefficient that will determine the persistence of the population assuming no restriction on the length of the patch. Here we use their estimated means and standard deviation of and to fix the values of , , , and defined above and
unless otherwise stated. We then study the persistence of the population in our stochastic framework for a fixed patch size km.
4.1 Critical shifting speed for Gaussian kernel
First, we investigate numerically the variation of the critical shifting speed as a function of the variance of the dispersal kernel in the case of Gaussian dispersal kernel. Indeed, using the analysis of Section 3.2, we know that for Gaussian kernel, when , that is the shifting speed is not stochastic anymore, then there exists a critical speed such that for all the population persists in the moving frame whereas if the population goes extinct. Using equation (3.46) we know that with being the principal eigenvalue of the linear operation defined in (3.31) when and being the geometric average of the geometric growth rate, defined in (3.28). Notice from the derivation of in the previous section that when