SIAM Journal on Applied Mathematics, 74 - University of Louisville ...
SIAM J. APPL. MATH. Vol. 74, No. 5, pp. 1397–1417
c 2014 Society for Industrial and Applied Mathematics
PERSISTENCE AND SPREAD OF A SPECIES WITH A SHIFTING HABITAT EDGE∗ BINGTUAN LI† , SHARON BEWICK‡ , JIN SHANG† , AND WILLIAM F. FAGAN‡ Abstract. We study a reaction-diffusion model that describes the growth and spread of a species along a shifting habitat gradient on which the species’ growth increases. It is assumed that the linearized species growth rate is positive near positive infinity and is negative near negative infinity. We show that the persistence and spreading dynamics depend on the speed of the shifting habitat edge c and a number c∗ (∞) that is determined by the maximum linearized growth rate and the diffusion coefficient. We demonstrate that if c > c∗ (∞), then the species will become extinct in the habitat, and that if c < c∗ (∞), then the species will persist and spread along the shifting habitat gradient at an asymptotic spreading speed c∗ (∞). Key words. reaction-diffusion equation, shifting habitat edge, persistence, spreading speed AMS subject classifications. 92D40, 92D25 DOI. 10.1137/130938463
1. Introduction. Whether a species can persist and spread in an environment is a fundamental and long-standing question in ecology. Recently, this question has been brought to the forefront by anthropogenic disturbance, including climate warming and landscape conversion. In particular, anthropogenic change has been blamed for widespread population declines across a broad range of taxa [30]. However, at the same time, there is evidence that other populations, most notably pest species and alien invasives, may be benefitting from recent changes [7, 24, 28]. Why certain species suffer while others tolerate and even prosper under changing conditions remains an area of active investigation. Early on, it was recognized that species persistence and spread often depends on spatial context, specifically the combined effects of habitat range and species dispersal. This led to a vast body of literature analyzing the role of spatial processes in population dynamics [12, 18]. Many of the conclusions from these early studies pertain to species persistence in changing landscapes as well. Indeed, because drivers like climate change alter the spatial spread of suitable habitats, population level effects can often be understood in terms of perturbations to spatial processes. While several frameworks exist for studying spatial population ecology, much analysis is based on reaction-diffusion models like the Fisher’s equation (1.1)
∂2u ∂u = d 2 + ru − u2 . ∂t ∂x
Here u is the population density of the species of interest, d is the diffusion coefficient, r is the population growth rate, and −u2 accounts for density-dependent death (e.g., ∗ Received by the editors September 25, 2013; accepted for publication (in revised form) May 29, 2014; published electronically September 16, 2014. http://www.siam.org/journals/siap/74-5/93846.html † Department of Mathematics, University of Louisville, Louisville, KY 40059 (bing.li@louisville. edu, [email protected]). The research of these authors was partially supported by the National Science Foundation under grant DMS-1225693. ‡ Department of Biology, University of Maryland, College Park, MD 20742 (sharon bewick@ hotmail.com, [email protected]). The research of these authors was partially supported by the National Science Foundation under grant DMS-1225917.
1397
1398
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
due to resource limitation). Frequently, models such as (1.1) are used to study the spreading speed of a population. Loosely speaking, this is the asymptotic rate at which a species, initially introduced in a bounded domain, expands its spatial range. Ecologically, an understanding of spreading speeds can provide insight into invasion processes, for example, how quickly an introduced species can move into a novel landscape or how rapidly an extirpated species can recover to its previous √ range. It is well known that for the Fisher’s equation, the spreading speed is c∗ = 2 dr, which is the slowest speed of a class of traveling wave solutions to (1.1) [12, 18]. A traveling wave solution describes the propagation of a species as a wave with a fixed shape and a fixed speed. A number of extensions of the basic Fisher’s equation exist. It has, for example, been studied in systems with general growth functions [2], in multiple dimensions [3], and in systems with interacting species [10, 19, 21, 48]. Modified versions have been applied to discrete time systems in order to capture seasonal life-histories. Finally, a number of different forms of spatial heterogeneity have been considered [1, 11, 39, 48, 50]. To incorporate spatial heterogeneity, one or both of the diffusion coefficient, d, or the population growth rate, r, are taken to be functions of space. Shigesada, Kawasaki, and Teramoto, for example, studied spreading speeds in systems where both the diffusion coefficient and the growth rate are periodic [39]. Weinberger considered a more general problem for both continuous time and discrete time models [49]. In both papers, periodicity was attributed to patchiness in the environment. Spatial heterogeneity has also been considered for models in time almost periodic and space periodic media [16]. Recently, there has been renewed interested in spatial ecology, largely driven by the threats associated with global change. In particular, due to the processes of global change, habitat regions suitable for population growth of many species have shifted geographically and are predicted to shift even more. For example, changes in the boundaries between forest and non-forest habitat types are receiving extensive study by ecologists and biogeographers in biomes worldwide (e.g., [14, 37, 38, 44]). In many cases, the boundary between suitable and unsuitable habitat shifts gradually over a period of years or decades (e.g., [31, 32]), but in other cases the range edge can shift quickly [37] or even in response to a single climatic event [40]. Changes in habitat boundaries may result in shrinkage of particular habitat types, such as when forest habitat expands upslope, encroaching on alpine communities [14, 25]. Similarly, expansion of scrub-woodland habitats into areas historically dominated by grasslands decreases the spatial extent of the grasslands [14, 31]. In other situations, suitable habitat is actually expanding [14, 25]. For example, along the Antarctic Peninsula, temperatures have been increasing rapidly for several decades, resulting in a decline in spring sea ice. This has increased accessibility of ice-free territory on shore, allowing the gentoo penguin (Pygoscelis papua) to spread southward beyond its historical limits [22]. Understanding and predicting the effects of changes like those described above has required the development of new theoretical models. In general, the emphasis of these models is on the role of changing habitat suitability. Thus while the reactiondiffusion framework in (1.1) is still valid, it requires population growth rates that are functions of time and location. One of the most common effects of global change is shifts or translations in habitat ranges. This is often modeled by assuming that the population growth rate can be described as r(x − ct), where r(x) is the spatially varying baseline, or historic rate of population growth. Initially suitable habitat is all locations with r(x) > 0. As global change progresses, however, the spatial locations of
SPREAD WITH A SHIFTING HABITAT EDGE
1399
the suitable habitats move in the positive direction at a rate of c units per timestep. Berestycki et al. [5] used a model like this to define the minimum size of a moving habitat necessary to sustain a non-zero species population (analogous to the critical patch size in models with stationary habitats). They then generalized their analysis to derive criteria for the persistence of a species in any region with a moving and spatially varying habitat. In a similar study, Potapov and Lewis [33] examined the implications of a moving habitat on competition between two species. More recently, Zhou and Kot [51] have extended the work in Berestycki et al. [5] to consider discrete time systems. Today, most applications of diffusion equations to global change scenarios have focused on determining criteria for the persistence of an established species. These are species that exist at an equilibrium distribution at the onset of climate change but then must track climate change in order to persist into the future. Fewer studies have considered the spread or invasion of an introduced species and how this might be impacted by shifts in habitat suitability. In [33], the authors consider the invasion problem from the perspective of two competing species and a moving boundary. Here, we take a broader and more basic perspective and consider the spread of a single species over a region with varying habitat suitability that is shifting in time. This scenario is relevant to a number of different problems in ecology, for example, predicting the spread of an invasive species under climate change [15, 41] or ascertaining the likelihood of success of assisted migration (i.e., deliberate translocation of species to areas with more suitable conditions) [23]. Indeed, understanding how climate change will affect non-equilibrial populations like those associated with invasion processes is a broad class of problems—one that has garnered some attention from biological communities [15, 34, 41] but has been largely ignored from a mathematical modeling perspective (but see Jeschke and Strayer [17]). To explore the issue of species spread in the context of climate change we consider the following formulation of the reaction-diffusion model in (1.1): (1.2)
∂2u ∂u = d 2 + ur(x − ct) − u2 . ∂t ∂x
Specifically, we assume an infinite domain −∞ < x < +∞ and a constant diffusion coefficient, d. For population growth, we take r(x − ct), where c > 0 and r(ξ) is continuous and nondecreasing and bounded with r(−∞) < 0 and r(∞) > 0. r(x − ct) thus divides the spatial domain into two parts: the region with good-quality habitat suitable for growth (i.e., r(x − ct) > 0), and the region with poor-quality habitat unsuitable for growth (i.e., r(x − ct) < 0). The edge of the habitat suitable for species growth is shifting at a speed c. We show that the persistence and spreading dynamics depend on c and c∗ (∞) = 2 dr(∞). We demonstrate that if c > c∗ (∞), then the species will become extinct in the habitat, and that if c < c∗ (∞), then the species will persist and spread along the shifting habitat gradient at an asymptotic spreading speed c∗ (∞). This paper is organized as follows. In the next section, the mathematical results regarding the spatial dynamics of (1.2) are presented. Section 3 presents a number of simulations that help to illustrate the results from section 2. Finally, in section 4 we discuss the relevance and implications of our results. 2. Mathematical results. In this section we provide mathematical results on the spatial dynamics for model (1.2). We begin with making the following hypothesis for r(x).
1400
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
Hypothesis 2.1. r(x) is continuous, nondecreasing and bounded, and piecewise continuously differentiable in x for −∞ < x < ∞, 0 < r(∞) < ∞, and −∞ < r(−∞) < 0. For 0 ≤ u1 , u2 ≤ r(∞), −∞ < x < ∞, and t ≥ 0, |u1 (r(x − ct) − u1 ) − u2 (r(x − ct) − u2 )| ≤ 3r(∞)|u1 − u2 |, so that u(r(x − ct) − u) is Lipschitz continuous in u. Clearly u ≡ 0 is a trivial solution (thus a lower solution) and u ≡ r(∞) is an upper solution of (1.2). The theory on the existence and uniqueness of solutions for reaction-diffusion systems has been well established (e.g., [27, Theorem 2.1], [46, Lemma 1.2]). It is known that the initial value problem of (1.2) with u(0, x) = u0 (x) where u0 (x) is continuous and 0 ≤ u0 (x) ≤ r(∞) has a unique classical solution u(t, x) with 0 ≤ u(t, x) ≤ r(∞). Let ρ > 3r(∞). Then u(ρ + r(x − ct) − u) is nondecreasing in u for 0 ≤ u ≤ r(∞). Equation (1.2) can be written as ∂2u ∂u + ρu = d 2 + u(ρ + r(x − ct) − u). ∂t ∂x
(2.1)
The solution of (2.1) with u(0, x) = u0 (x) satisfies the integral equation +∞ u(t, x) = k(t, x − y)u0 (y)dy −∞ t
+∞
+ 0
−∞
k(t − τ, x − y)u(τ, y)[ρ + r(y − cτ ) − u(τ, y)]dydτ,
where y2 1 e−ρs− 4ds . k(s, y) = √ 4πds
Consider the sequence u(n) (t, x) defined by (2.2) u
(n+1)
(t, x) =
+∞
−∞ t
k(t, x − y)u0 (y)dy
+∞
+ 0
−∞
k(t − τ, x − y)u(n) (τ, y)[ρ + r(y − cτ ) − u(n) (τ, y)]dydτ,
where u(0) (t, x) = 0 or u(0) (t, x) = r(∞). Lemma 7.22 in [27] shows that if u(0) (t, x) = 0, then u(n) (t, x) is nondecreasing in n, and if u(0) (t, x) = r(∞), then u(n) (t, x) is nonincreasing in n, 0 ≤ u(n) (t, x) ≤ r(∞). In both cases, u(t, x) = limn→∞ u(n) (t, x) is the solution of (2.1) with u(0, x) = u0 (x), and 0 ≤ u(t, x) ≤ r(∞). We first provide a useful lemma for the equation (2.3)
∂2u ∂u = d 2 + u(r(x) − u). ∂t ∂x
Lemma 2.1. Let u ¯(t, x) be the solution of (2.3) with u ¯(0, x) = r(∞). Then u ¯(t, x) ¯(t, −∞) = 0, and u ¯(t, ∞) ≡ is nonincreasing in t and nondecreasing in x, limt→∞ u r(∞) for t > 0.
SPREAD WITH A SHIFTING HABITAT EDGE
1401
Proof. It follows from (2.2) with u(0) (t, x) = r(∞) and c = 0 that the solution u ¯(t, x) of (2.3) with u ¯(0, x) = r(∞) has the property that u ¯(t, x) = limn→∞ u(n) (t, x), where (2.4) u(n+1) (t, x) = e−ρt r(∞) t +∞ + k(s, y)u(n) (t − s, x − y)[ρ + r(x − y) − u(n) (t − s, y)]dyds. 0
−∞
Here we have used the simple fact that +∞ k(t, y)dy = e−ρt . −∞
We first consider n = 1 in (2.4). Direct calculations show that +∞ ∂u(1) (t, x) k(t, y)r(∞)[ρ + r(x − y) − r(∞)]dy = −ρe−ρt r(∞) + ∂t −∞ +∞ k(t, y)r(∞)[r(x − y) − r(∞)]dy = −∞
≤0 and ∂u(1) (t, x) = ∂x Assume that
∂u(n) (t,x) ∂t
t 0
≤ 0 and
+∞
−∞
k(s, y)r(∞)r (x − y)dyds ≥ 0.
∂u(n) (t,x) ∂x
≥ 0 for some n > 0. Then
∂u(n+1) (t, x) = −ρe−ρt r(∞) ∂t +∞ k(t, y)u(n) (0, x − y)[ρ + r(x − y) − u(n) (0, x − y)]dy + −∞
t
+∞
+ 0
−∞
k(s, y)
∂u(n) (t, x) [ρ + r(x − y) − 2u(n) (t − s, y)]dyds. ∂t
Note that ρ + r(x − y) − 2u(n) (t − s, y) ≥ 0 because of ρ > 3r(∞) and that 0 ≤ u(n) (t − s, y) ≤ r(∞). Since (2.4) shows u(n) (0, x) = r(∞) for all n, we have that +∞ ∂u(n+1) (t, x) = k(t, y)r(∞)[r(x − y) − r(∞)]dy ∂t −∞ t +∞ ∂u(n) (t, x) [ρ + r(x − y) − 2u(n) (t − s, y)]dyds k(s, y) + ∂t 0 −∞ ≤ 0. On the other hand, (n) t +∞ ∂u (t − s, x − y) ∂u(n+1) (t, x) = (ρ + r(x − y) − 2u(n) (t − s, x − y)) k(s, y) ∂x ∂x 0 −∞ (n) + r (x − y)u (t − s, x − y) dyds ≥ 0.
1402
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN (n)
By induction, ∂u limn→∞ u(n) (t, x),
(n)
(t,x) ∂t ∂u ¯ (t,x) ∂t
(2.5) u ¯(t, x) = r(∞)e−ρt +
≤ 0 and ∂u ∂x(t,x) ≥ 0 for all n > 0. Since u ¯(t, x) = ∂u ¯ (t,x) ≤ 0 and ∂x ≥ 0. u ¯ satisfies
t 0
+∞
−∞
k(s, y)¯ u(t − s, x − y)[ρ + r(x − y) − u ¯(t − s, x − y)]dyds.
Taking the limit x → −∞ in (2.5) and applying the dominated convergence theorem we obtain t +∞ ¯(t − s, −∞)] k(s, y)dyds u ¯(t, −∞) = r(∞)e−ρt + u¯(t − s, −∞)[ρ + r(−∞) − u 0
= r(∞)e−ρt + = r(∞)e−ρt +
−∞
t
0
0
t
e−ρs u ¯(t − s, −∞)[ρ + r(−∞) − u ¯(t − s, −∞)]ds e−ρ(t−s) u ¯(s, −∞)[ρ + r(−∞) − u ¯(s, −∞)]ds.
This implies that u ¯(t, −∞) satisfies the differential equation ∂u ¯(t, −∞) =u ¯(t, −∞)[r(−∞) − u ¯(t, −∞)]. ∂t Since r(−∞) < 0, 0 is the only nonnegative equilibrium of the differential equation and it attracts all nonnegative solutions. We therefore have limt→∞ u ¯(t, −∞) = 0. Similarly one can show that u ¯(t, ∞) satisfies the differential equation ∂u ¯(t, ∞) =u ¯(t, ∞)[r(∞) − u ¯(t, ∞)], ∂t which has two equilibria 0 and r(∞). Since u ¯(0, ∞) = r(∞), we have that u ¯(t, ∞) ≡ r(∞). The proof is complete. For r(x) > 0, define c∗ (x) = 2 dr(x). It is easily seen that c∗ (x) = inf φ(x; μ), μ>0
where φ(x; μ) = The infimum occurs at μ∗ (x) =
r(x) d .
dμ2 + r(x) . μ
The function
ψ(μ) = 2dμ is useful. It is easily seen that φ(x; μ) > ψ(μ) for 0 < μ < μ∗ (x) and φ(x; μ∗ (x)) = ψ(μ∗ (x)).
SPREAD WITH A SHIFTING HABITAT EDGE
1403
Our first theorem shows that if c∗ (∞) < c, then the species will eventually become extinct in space under certain conditions. Theorem 2.1. Assume that Hypothesis 2.1 is satisfied. Let c > 0 and c∗ (∞) < c. If 0 ≤ u0 (x) ≤ r(∞) and u0 (x) ≡ 0 for all sufficiently large x, then for every ε > 0 there exits T > 0 such that for t ≥ T , the solution u(t, x) of (1.2) with u(0, x) = u0 (x) satisfies u(t, x) ≤ ε for all x. Proof. Let u ¯(t, x) be the solution of (2.3) with u¯(0, x) = r(∞). Since ¯(t, −∞) = 0 by Lemma 2.1, for every ε > 0, there exist T1 > 0 and M > 0 limt→∞ u such that for x ≤ −M , u¯(T1 , x) < ε. Since u ¯(t, x) is nonincreasing in t, (2.6)
u ¯(t, x) < ε for t ≥ T1 and x ≤ −M.
The function u ˜(t, x) = u ¯(t, x − ct) satisfies (2.7)
∂u ˜ ∂2u ˜ ∂u ˜ =d 2 +c +u ˜(r(x − ct) − u ˜). ∂t ∂x ∂x
Since u ¯(t, x) is nondecreasing in x,
∂u ˜ ∂x
≥ 0. Since c > c∗ (∞) > 0, (2.7) shows that
∂ 2 u˜ ∂u ˜ ≥ d 2 + u˜(r(x − ct) − u ˜), ∂t ∂x so that u˜(t, x) = u¯(t, x − ct) with u˜(0, x) = r(∞) is an upper solution of (1.2). Since u(0, x) = u0 (x) ≤ r(∞), u(t, x) ≤ u ¯(t, x − ct). It follows from this and (2.6) that (2.8)
u(t, x) < ε for t ≥ T1 and x ≤ −M + ct.
Choose 0 < δ < c − c∗ (∞). Let μδ > 0 be the smaller positive solution of ∗ φ(∞; μ) = c∗ + δ/2. It is easily seen that u ˆ(t, x) = Ae−μδ (x−(c (∞)+δ/2)t) , where A is a positive constant, is a solution of the linear equation ∂2u ∂u = d 2 + r(∞)u. ∂t ∂x Since r(∞)u ≥ u(r(x − ct) − u), u ˆ(t, x) is an upper solution of (1.2). Choose A sufficiently large such that u0 (x) ≤ uˆ(0, x) = Ae−μδ x . Then u(t, x) ≤ Ae−μδ (x−(c
∗
(∞)+δ/2)t)
.
This shows that for x ≥ (c∗ (∞) + δ)t, u(t, x) ≤ Ae−(μδ δ/2)t . It follows that for the above given ε > 0, there exists T2 > 0 such that (2.9)
u(t, x) < ε for t ≥ T2 and x ≥ (c∗ (∞) + δ)t.
On the other hand, since c > c∗ (∞) + δ, there exists T3 > 0 such that for t ≥ T3 , −M +ct > (c∗ (∞)+δ)t. This, (2.8), and (2.9) show that for t ≥ T := max{T1 , T2 , T3 } and for all x, u(t, x) < ε. The proof is complete. Definition 2.2. We call a function u a continuous weak lower solution of model (1.2) if u is continuous for t ≥ 0 and −∞ < x < ∞, and ∂2u ∂u ≤ d 2 + u(r(x − ct) − u) ∂t ∂x
1404
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
in the distributional sense, i.e., for any T > 0 and any η ∈ C 2,1 ((−∞, ∞) × [0, T ]) with η ≥ 0 and supp η(·, t) being a bounded interval for all t ∈ [0, T ],
∞
−∞
≤
1 u(t, x)η(t, x)dx|t=T t=0
T1
∞
−∞
0
[u(s, x)(dηxx + ηt )(s, x) + η(s, x)u(s, x)(r(x − cs) − u(s, x))]dxds
if T1 ∈ [0, T ]. This definition is a slightly modified version of Definition 1.1 in [46]. Weak lower solutions were used in [2, 46] in studying reaction-diffusion systems. For γ > 0 and μ > 0, define (2.10)
v(μ; x) =
e−μx sin γx
if 0 ≤ x ≤ π/γ,
0
elsewhere.
(A scalar multiple of this function is denoted by v(s) in Weinberger [47].) v(μ; x) is a continuous function in x and its second order derivative in x exists and is continuous when x = 0, π/γ. The maximum of v(μ; x) occurs at σ(μ) = (1/γ) tan−1 (γ/μ). σ(μ) is a strictly decreasing function of μ. We have the following useful lemma. ∗
Lemma 2.3. Assume that c∗ (∞) > c ≥ 0. For any satisfying 0 < < c (∞)−c , 3 let be a number such that c∗ () = c∗ (∞) − . Let 0 < μ1 < μ2 < μ∗ () with ψ(μ1 ) = c + and ψ(μ2 ) = c∗ (∞) − 2. Then for any μ ∈ [μ1 , μ2 ], there exist a > 0 and γ > 0 sufficiently small such that av(μ; x − − ψ(μ)t) with v given by (2.10) is a continuous weak lower solution of (1.2). Furthermore, if u(0, x) ≥ av(μ; x − ), then u(t, x) ≥ av(μ; x − − ψ(μ)t) for all t > 0. Proof. By using the definition of v and integration by parts, we find that
T1
+∞
(2.11)
v(μ; x − − zs)(dηxx + ηt )(s, x)dxds
−∞ T1 +zs+π/γ
0
= 0
+zs T1 +zs+π/γ
0
+zs +zs+π/γ
=
+ +zs T1
+γ
0
e
v(μ; x − − zs)(dηxx + ηt )(s, x)dxds [dvxx (μ; x − − zs) − vs (μ; x − − zs)]η(s, x)dxds
v(μ; x − − zs)η(s, x)|T0 1 dxds
−μ π γ
π η s, + zs + + η(s, + zs) ds. γ
Direct calculations show that for μ ∈ [μ1 , μ2 ], x = + ψ(μ)t, and x = + ψ(μ)t + π/γ, (2.12)
dvxx (μ; x − − ψ(μ)t) − vt (μ; x − − ψ(μ)t) = −d(μ2 + γ 2 )v(μ; x − − ψ(μ)t).
SPREAD WITH A SHIFTING HABITAT EDGE
1405
Let f (t, x, u) = u(r(x − ct) − u). It follows from (2.11) and (2.12) that for μ ∈ [μ1 , μ2 ] and for sufficiently small a and γ, u ˜(μ; t, x) = av(μ; x − − ψ(μ)t) satisfies 0
T1
+∞
[˜ u(μ; s, x)(dηxx + ηt )(s, x) + η(s, x)f (s, x, u˜(μ; s, x))]dxds
−∞ ∞
−
=a
1 u˜(μ; s, x)η(s, x)dx|t=T t=0
−∞ T1 +ψ(μ)s+π/γ
0
[r(x − cs) − dμ2 − dγ 2 − av(μ; x − − ψ(μ)s)]
+ψ(μ)s
× v(μ; x − − ψ(μ)s)η(s, x)dxds T1 π −μ π γ + aγ e η s, + zs + + η(s, + zs) ds γ 0 T1 +ψ(μ)s+π/γ [r() − dμ22 − dγ 2 − av(μ; x − − ψ(μ)s)] ≥a 0
+ψ(μ)s
× v(μ; x − − ψ(μ)s)η(s, x)dxds T1 π −μ π γ η s, + zs + + aγ e + η(s, + zs) ds γ 0 T1 +ψ(μ)s+π/γ
∗ (c (∞) − 1.5) − dγ 2 − av(μ; x − − ψ(μ)s) =a 2d +ψ(μ)s 0 × v(μ; x − − ψ(μ)s)η(s, x)dxds T1 π −μ π γ + aγ η s, + zs + e + η(s, + zs) ds γ 0 T1 +ψ(μ)s+π/γ (c + ) − dγ 2 − av(μ; x − − ψ(μ)s) ≥a 2d +ψ(μ)s 0 × v(μ; x − − ψ(μ)s)η(s, x)dxds T1 π −μ π γ η s, + zs + + aγ e + η(s, + zs) ds γ 0 ≥ 0. It follows from Definition 2.2 that for μ ∈ [μ1 , μ2 ] and for sufficiently small a > 0 and γ > 0, av(μ; x − − ψ(μ)t) is a continuous weak lower solution of (1.2). If u(0, x) ≥ av(μ; x − ), then it follows from Lemma 1.2 in Wang [46] that u(t, x) ≥ av(μ; x − − ψ(μ)t) for all t > 0. The proof is complete. The following theorem shows that if c∗ (∞) > c, then the species persists in space and spreads to the right at the asymptotic spreading speed c∗ (∞). Theorem 2.2. Assume that Hypothesis 2.1 is satisfied. Let c∗ (∞) > c ≥ 0. Then the following statements are valid: (i) If 0 ≤ u(0, x) ≤ r(∞), then for any ε > 0, lim sup u(t, x) = 0. t→∞
x≤t(c−ε)
(ii) If 0 ≤ u(0, x) ≤ r(∞), and u(0, x) ≡ 0 for all sufficiently large x, then for any ε > 0,
1406
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
lim
t→∞
sup x≥t(c∗ (∞)+ε)
u(t, x) = 0.
(iii) If 0 ≤ u(0, x) ≤ r(∞), and u(0, x) > 0 on a closed interval, then for every ε with 0 < ε < (c∗ (∞) − c)/2, |r(∞) − u(t, x)| = 0. sup lim t→∞
t(c+ε)≤x≤t(c∗ (∞)−ε)
Proof. The first part of the proof of Theorem 2.1 shows that u ¯(t, x − ct) with u ¯(0, x) = r(∞) is an upper solution of (1.2). Since u(0, x) ≤ r(∞), u(t, x) ≤ u ¯(t, x − ct). Lemma 2.1 shows that for any > 0, there exist T1 > 0 and a real number M such that for t ≥ T1 and x ≤ −M + ct, u(t, x) < . On the other hand, for any ε > 0 there exists T2 > 0 such that for t ≥ T2 , (c − ε)t ≤ −M + ct. We therefore have that for t ≥ max{T1 , T2 } and x ≤ t(c − ε), u(t, x) < . This proves statement (i). The proof of statement (ii) is similar to the second half of the proof of Theorem 2.1 and is omitted. We now prove statement (iii). For any small positive with 0 < < ∗ }, let be a number such that min{ 13 , r(∞), c (∞)−c 3 c∗ () = c∗ (∞) − . Choose 0 < μ1 < μ2 < μ∗ () such that ψ(μ1 ) = c + and ψ(μ2 ) = c∗ (∞) − 2. By α v(μ; x−−ψ(μ)t) Lemma 2.3, for any μ ∈ [μ1 , μ2 ] and small α > 0 and γ > 0, v(μ;σ(μ)) with v given by (2.10) is a continuous weak lower solution of (1.2). Since u(0, x) ≥ 0 and u(0, x) ≡ 0, u(t, x) > 0 for all x and any t > 0. Choose σ(μ1 ) 0 < t0 < ψ(μ , and choose α and γ sufficiently small such that u(t0 , x) ≥ α for 1) x ∈ [, + 4π/γ]. Define ⎧ α ⎪ v(μ1 ; x − ) if ≤ x ≤ + σ(μ1 ), ⎪ ⎪ ⎪ v(μ ; 1 σ(μ1 )) ⎪ ⎪ ⎪ ⎨ α if + σ(μ1 ) ≤ x ≤ + 3π/γ + σ(μ2 ), w(0, x) = α ⎪ ⎪ v(μ2 ; x − − 3π/γ) if + 3π/γ + σ(μ2 ) ≤ x ≤ + 4π/γ, ⎪ ⎪ v(μ2 ; σ(μ2 )) ⎪ ⎪ ⎪ ⎩ 0 elsewhere. α v(μ1 ; x − − s) for 0 ≤ s ≤ 2π/γ, and It is easily seen that w(0, x) ≥ v(μ1 ;σ(μ 1 )) α w(0, x) ≥ v(μ2 ;σ(μ2 )) v(μ2 ; x − − 3π/γ + s) for 0 ≤ s ≤ 2π/γ. Since u(t0 , x) ≥ α for x ∈ [, + 4π/γ], Lemma 2.2 shows that for t ≥ t0 and 0 ≤ s ≤ 2π/γ,
(2.13)
u(t, x) ≥
α v(μ1 ; x − − ψ(μ1 )(t − t0 ) − s) v(μ1 ; σ(μ1 ))
SPREAD WITH A SHIFTING HABITAT EDGE
1407
and (2.14)
u(t, x) ≥
α v(μ2 ; x − − 3π/γ − ψ(μ2 )(t − t0 ) + s). v(μ2 ; σ(μ2 ))
Equation (2.13) implies that for t ≥ t0 , (2.15) u(t, x) ⎧ α ⎪ if + ψ(μ1 )(t − t0 ) ≤ x ≤ + σ(μ1 ) ⎪ ⎪ v(μ1 ; σ(μ1 )) ⎪ ⎪ ⎪ ⎪ ⎪ + ψ(μ1 )(t − t0 ), ⎨ × v(μ1 ; x − − ψ(μ1 )(t − t0 )) ≥ α if + σ(μ1 ) + ψ(μ1 )(t − t0 ) ≤ x ≤ + 2π/γ ⎪ ⎪ ⎪ ⎪ ⎪ + σ(μ1 ) + ψ(μ1 )(t − t0 ), ⎪ ⎪ ⎪ ⎩ 0 elsewhere. On the other hand, (2.14) indicates that for t ≥ t0 , (2.16) u(t, x) ⎧ α ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ α ≥ v(μ ; ⎪ 2 σ(μ2 )) ⎪ ⎪ ⎪ ⎪ × v(μ2 ; x − − 3π/γ − ψ(μ2 )(t − t0 )) ⎪ ⎪ ⎪ ⎩ 0
if + π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) ≤ x ≤ + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ), if + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) ≤ x ≤ + 4π/γ + ψ(μ2 )(t − t0 ), elsewhere.
1 )−σ(μ2 ) Let h = π/γ+σ(μ ψ(μ2 )−ψ(μ1 ) . Since + 2π/γ + σ(μ1 ) + ψ(μ1 )(t − t0 ) ≥ + π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) for t0 ≤ t ≤ t0 + h, (2.15) and (2.16) show that for t0 ≤ t ≤ t0 + h,
(2.17)
u(t, x) ≥ w(t − t0 , x),
where (2.18) w(t − t0 , x) ⎧ α ⎪ ⎪ ⎪ ⎪ v(μ ; 1 σ(μ1 )) ⎪ ⎪ ⎪ ⎪ ⎪ × v(μ1 ; x − − ψ(μ1 )(t − t0 )) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ α = ⎪ ⎪ α ⎪ ⎪ ⎪ ⎪ ⎪ v(μ ; 2 σ(μ2 )) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ × v(μ2 ; x − − 3π/γ − ψ(μ2 )(t − t0 )) ⎪ ⎪ ⎩ 0
if + ψ(μ1 )(t − t0 ) ≤ x ≤ + σ(μ1 ) + ψ(μ1 )(t − t0 ), if + σ(μ1 ) + ψ(μ1 )(t − t0 ) ≤ x ≤ + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ), if + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) ≤ x ≤ + 4π/γ + ψ(μ2 )(t − t0 ), elsewhere.
1408
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
We claim that (2.17) is valid for all t ≥ t0 . Assume that (2.17) is true for t0 ≤ α t ≤ t0 + nh for some positive integer n. Then w(nh, x) ≥ v(μ1 ;σ(μ v(μ1 ; x − − 1 )) α v(μ ; x − − 3π/γ − nhψ(μ ) nhψ(μ1 ) − s) and w(nh, x) ≥ v(μ2 ;σ(μ 2 2 + s) for 0 ≤ 2 )) s ≤ 2π/γ + (ψ(μ2 ) − ψ(μ1 ))nh. It follows from Lemma 2.3 that for t ≥ t0 + nh and 0 ≤ s ≤ 2π/γ + (ψ(μ2 ) − ψ(μ1 ))nh, (2.19)
u(t, x) ≥
α v(μ1 ; x − − nhψ(μ1 ) − ψ(μ1 )(t − (t0 + nh)) − s) v(μ1 ; σ(μ1 ))
and (2.20) u(t, x) ≥
α v(μ2 ; x − − 3π/γ − nhψ(μ2 ) − ψ(μ2 )(t − (t0 + nh)) + s). v(μ2 ; σ(μ2 ))
Equation (2.19) shows that u(t, x) ≥ α
(2.21) for
t ≥ t0 + nh
(2.22) and (2.23)
+ σ(μ1 ) + nhψ(μ1 ) + ψ(μ1 )(t − (t0 + nh)) ≤ x ≤ + 2π/γ + σ(μ1 ) + nhψ(μ2 ) + ψ(μ1 )(t − (t0 + nh)).
Equation (2.20) implies that (2.21) holds if (2.22) is satisfied and if (2.24)
+ π/γ + σ(μ2 ) + nhψ(μ1 ) + ψ(μ2 )(t − (t0 + nh)) ≤ x ≤ + 3π/γ + σ(μ2 ) + nhψ(μ2 ) + ψ(μ2 )(t − (t0 + nh)).
The intervals described by (2.23) and (2.24) overlap if + 2π/γ + σ(μ1 ) + nhψ(μ2 ) + ψ(μ1 )(t − (t0 + nh)) ≥ + π/γ + σ(μ2 ) + nhψ(μ1 ) + ψ(μ2 )(t − (t0 + nh)). This is equivalent to t ≤ t0 + 2nh +
π/γ + σ(μ1 ) − σ(μ2 ) = t0 + (2n + 1)h. ψ(μ2 ) − ψ(μ1 )
It follows that the intervals described by (2.23) and (2.24) overlap if t0 + nh ≤ t ≤ t0 + (n + 1)h. We therefore have that (2.21) holds if t0 + nh ≤ t ≤ t0 + (n + 1)h and + σ(μ1 ) + nhψ(μ1 ) + ψ(μ1 )(t − (t0 + nh)) ≤ x ≤ + 3π/γ + σ(μ2 ) + nhψ(μ2 ) + ψ(μ2 )(t − (t0 + nh)). It follows from this, (2.19), and (2.20) that (2.17) holds for t0 + nh ≤ t ≤ t0 + (n + 1)h. By induction, (2.17) holds for all t ≥ t0 . For the chosen > 0, there is L > 0 such that L 2 1 √ e−x dx ≥ 1 − . π −L
SPREAD WITH A SHIFTING HABITAT EDGE
1409
Note that for any s > 0, √ L 4ds
√ −L 4ds
x2 1 √ e− 4ds dx = 4πds
L
−L
2 1 √ e−ξ dξ. π
Let t1 > t0 be a sufficiently large number. For t > t1 , the solution u(t, x) satisfies the integral equation +∞ (2.25) k(t − t1 , x − y)u(t1 , y)dy u(t, x) = −∞ t
+∞
+ t1
−∞
k(t − τ, x − y)u(τ, y)[ρ + r(y − cτ ) − u(τ, y)]dydτ.
It follows from this and (2.17) that for t > t1 , (2.26)
u(t, x) ≥
+∞
−∞ t
k(t − t1 , x − y)w(t1 − t0 , y)dy
+∞
+ −∞
t1
k(t − τ, x − y)w(τ − t0 , y)[ρ + r(y − cτ ) − w(τ − t0 , y)]dydτ.
For t > t1 , x satisfying (2.27)
+ σ(μ1 ) + ψ(μ1 )(t − t0 ) + L 4d(t − t1 ) ≤ x ≤ + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) − L
and y satisfying
4d(t − t1 ),
−L 4d(t − t1 ) ≤ y ≤ L 4d(t − t1 ),
(2.28) we have that
+ σ(μ1 ) + ψ(μ1 )(t − t0 ) ≤ x − y ≤ + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) and x − y − ct ≥ + σ(μ1 ) + ψ(μ1 )(t − t0 ) − ct = + t + σ(μ1 ) − ψ(μ1 )t0 > . It follows from this and (2.18) that for x satisfying (2.27), +∞ (2.29) k(t − t1 , x − y)w(t1 − t0 , y)dy −∞
+∞
k(t − t1 , y)w(t1 − t0 , x − y)dy L√4d(t−t1 ) y2 1 − −ρ(t−t1 ) ≥e e 4d(t−t1 ) w(t1 − t0 , x − y)dy √ 4πd(t − t1 ) −L 4d(t−t1 ) √ L 4d(t−t1 ) y2 1 − e 4d(t−t1 ) dy ≥ αe−ρ(t−t1 ) √ 4πd(t − t1 ) −L 4d(t−t1 ) =
−∞
≥ (1 − )αe−ρ(t−t1 )
1410
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
and
+∞
(2.30) −∞
k(t − τ, x − y)w(τ − t0 , y)[ρ + r(y − cτ ) − w(τ − t0 , y)]dy
≥ (1 − )α[ρ + r(∞) − − α]e−ρ(t−τ ) . Here we have used the fact that for x satisfying (2.27) and y satisfying (2.28) r(x − y − ct) ≥ r() > r(∞) − . It follows from (2.26), (2.29), and (2.30) that for t ≥ t1 and x satisfying (2.27) u(t, x) ≥ u ˜(1) (t), where u ˜
(1)
(t) = (1 − )αe
−ρ(t−t1 )
+ (1 − )
t
e−ρ(t−τ ) α[ρ + r(∞) − − α]dτ.
t1
Equation (2.25) and induction show that for t ≥ t1 and x satisfying + σ(μ1 ) + ψ(μ1 )(t − t0 ) + nL 4d(t − t1 ) (2.31) ≤ x ≤ + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) − nL 4d(t − t1 ), u(t, x) ≥ u ˜(n) (t),
(2.32) where u ˜(n) (t) satisfies (2.33)
u ˜(n) (t) = (1 − )αe−ρ(t−t1 ) t + (1 − ) e−ρ(t−τ ) u ˜(n−1) (τ )[ρ + r(∞) − − u˜(n−1) (τ )]dτ. t1
Direct calculations and induction show that (2.34)
u˜(n) (t) = an + bn (t)e−ρ(t−t1 ) ,
where an = (1 − )an−1 (ρ + r(∞) − − an−1 )/ρ, a1 = (1 − )α[ρ + r(∞) − − α]/ρ, and bn (t) is a sum of polynomials, and products of polynomials and exponential functions ˜(n) (t) = an . in the form of e−jρ(t−t1 ) with j a positive integer. Observe that limt→∞ u For small and α, an increases to r(∞) − − ρ/(1 − ) as n → ∞. This and (2.34) show that there exist a positive integer N and t2 > t1 such that for t > t2 and x satisfies (2.31) with n replaced by N , (2.35)
u ˜(N ) (t) ≥ r(∞) − 2 − ρ/(1 − ).
We choose t1 sufficiently large such that for t ≥ t1 , + σ(μ1 ) + ψ(μ1 )(t − t0 ) + N L 4d(t − t1 )
< + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) − N L 4d(t − t1 ).
For any given ε with 0 < ε < (c∗ (∞) − c)/2, choose sufficiently small such that < ε/2. Then there exists t3 > t2 such that for t > t3 , + σ(μ1 ) + ψ(μ1 )(t − t0 ) + N L 4d(t − t1 ) < t(c + ε) < t(c∗ (∞) − ε)
1411
SPREAD WITH A SHIFTING HABITAT EDGE
and t(c∗ (∞) − ε) < + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) − N L 4d(t − t1 ). This shows that for t > t3 , t(c + ε) ≤ x ≤ t(c∗ (∞) − ε) implies that x satisfies (2.31) with n replaced by N . It follows from (2.32) and (2.35) that u(t, x) ≥ r(∞) − 2 − ρ/(1 − ). (2.36) lim inf ∗ t→∞ t(c+ε)≤x≤t(c (∞)−ε)
Since is arbitrarily small and u(t, x) ≤ r(∞) for all x and t, (2.36) shows that statement (iii) holds. The proof is complete. 3. Simulations. In this section we present some numerical simulations to the model (1.2) with −0.5 if x ≤ ct, r(x − ct) = 1 elsewhere, where c > 0, and the initial data 0.5 sin(x − 20) if 20 ≤ x ≤ 20 + π, (3.1) u(0, x) = 0 elsewhere. r(x − ct) moves to the right at a speed c, and the species is initially introduced on the interval [20, 20 + π]. Choose d = 1 so that c∗ (∞) = 2. Numerical simulations were conducted using MATLAB. Figure 1 displays the numerical solution with c = 2.1 > c∗ (∞), which shows that the solution eventually becomes extinct. Figure 2 displays
1.2 t=35
Population density
1
t=90 0.8
0.6
t=1
0.4
0.2 t=165 0
0
50
100
150
200 250 Location(x)
300
350
400
Fig. 1. A numerical approximation to the graph of u(t, x) with c = 2.1 and the initial data given by (3.1). The solution finally becomes extinct.
1412
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
1.2 t=35
Population density
1
t=90
t=165
0.8
0.6 t=1 0.4
0.2
0
0
50
100
150
200 250 Location(x)
300
350
400
Fig. 2. A numerical approximation to the graph of u(t, x) with c = 1.5 and the initial data given by (3.1). The solution persists and spreads to the right at the asymptotic speed c∗ (∞) = 2.
the numerical solution with c = 1.5 < c∗ (∞), which shows that the solution expends its spatial range to the right at the asymptotic speed c∗ (∞) = 2. The numerical results support Theorems 2.1 and 2.2. 4. Discussion. The ability to persist in a changing landscape is one of the biggest challenges faced by species today. Indeed, because of global change, species are being forced to shift their historical ranges in order to keep up with shifts in suitable habitat. Most obviously, this involves movement poleward or to higher elevations. However, it can also include movement along moisture or precipitation gradients, salinity gradients, or other abiotic gradients impacted by global change. While the issue of persistence in a changing landscape has been addressed by a number of theoretical studies, these have primarily focused on established species that exist at equilibrium distributions in a bounded domain prior to the onset of climate change. In this case, species are restricted to and distributed over a band of suitable habitat. Persistence then requires that a species keep pace with the movement of its habitat band. This is a good model for species that are limited by abiotic constraints at both an upper and a lower boundary, as would be the case for native species restricted to a range of elevations or latitudes. A different scenario emerges for species that are established over a finite range but, at least for the forseeable future, only limited at one boundary. For this scenario, persistence requires that the species spread away from the boundary fast enough to outrun boundary encroachment. In contrast to the persistence of established species, this problem has received less attention. It is, however, a good model for invasive species that are initially introduced to a small portion of their suitable habitat. It is also a good model for native species that, because of some historical event, have not reached an equilibrium distribution across their native range. Examples include
SPREAD WITH A SHIFTING HABITAT EDGE
1413
certain tree species [42, 43], as well as a number of amphibian and reptile species [4] whose current distributions still show evidence of the Last Glacial Maximum (LGM). Finally, understanding when a species can outrun a moving boundary has implications for assisted migration [23]. If, for example, a species is moved from a shrinking habitat at the top of a mountain to a much larger habitat at a higher latitude, it would be useful to know whether the introduced population will persist, given the encroachment of unsuitable conditions from the south. Very generally, we find that a spreading species population will go extinct if its habitat boundary moves at speeds greater than the maximum rate of expansion of the species population, c > c∗ (∞). On the other hand, the population will persist if its habitat boundary moves at lower speeds c < c∗ (∞). Ultimately, this defines whether a species will be able to outrun an encroaching boundary. First, let us consider our results in the context of invasive species. A large number of invasives are currently cold-limited and thus are expected to benefit from climate change [7, 24, 28]. However, others will be negatively impacted. These species will lose habitat as a result of climate change, which may, in turn, limit their ability to persist and spread as alien invasives. This has led some authors to speculate that climate change will be a zero-sum game for invasive species [15]. In particular, while climate change will likely lead to the worsening or new emergence of a handful of invasives, it should reduce the impact of others, ultimately resulting in no change to the net burden of invasive species. Not surprisingly, less is known about species that show range contraction or reduced impacts as a result of climate change [45]. Nevertheless, several examples exist, for instance, the two fish species, sunbleak and ide, in England [8] as well as a number of Eurasion annual grasses in Southern Africa [29] and perennial grasses in Australia [13]. A key finding of our model is that despite losing habitat at one boundary, introduced species may continue spreading, provided that they can expand at a rate faster than climate change. In our model, the rate of expansion of a particular population depends on how well the species disperses, d, as well as its linearized maximum rate of population growth, r(∞). Interestingly, the majority of invasive species are characterized by rapid population growth and enhanced dispersal. In tree species, for instance, invasiveness has been associated with short juvenile periods, short intervals between large seed crops, and small seed size [35]. The first two traits lead to increased rates of population growth, while small seed size can increase dispersal. Consequently, most invasive species will be characterized by both large d and large r(∞), suggesting that climate change will rarely, if ever, halt the spread of an invasive species. This appears to be the case for some of the annual grasses in Southern Africa that are continuing their spread to high alpine communities with devastating effects [29]. The situation is different for native species with nonequilibrium distributions. This can be seen most dramatically with the example of species still spreading northward from ice-age refugia. Quite obviously, these species are expanding at a minimal rate. If this is because they are intrinsically poor at breeding and dispersal, then given estimates of current climate velocity on the order of 80 to 2700 m/yr [9, 20], it seems likely that species still strongly associated with ice-age refugia will have c > c∗ (∞). Consequently, existing populations of these species will not expand fast enough to escape encroachment by warmer climates. For northern hemisphere species with a southern range limit, the end result will be extinction. While a number of European species limited to ice-age refugia are Mediterranean and thus less likely to be severely limited by warmer climates, others are associated with alpine or Pannonian climate zones [26]. Species in these latter categories, for example, Pritzelago alpina or Pulsatilla alpina, are thought to have undergone range shifts and range
1414
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
contractions at their southern and/or low elevation boundaries in response to postglacial warming following the LGM [6, 26]. Consequently, these species are also likely to exhibit range losses in response to current warming regimes. Because they are limited by climate at their southern (low elevation) range boundary and by spreading speeds at their northern (high elevation) range boundary, we predict that alpine and Pannionian species that currently exhibit a strong association with ice-age refugia will be at a particularly high risk for extinction under ongoing global change. Moreover, transporting these species to higher latitudes will be of little use since, even there, they will ultimately be overtaken by climates that are unsuitable for growth. While the assisted migration of species from ice-age refugia may be futile, our model does not imply that all assisted migration will fail. Species with high dispersal and reproductive rates should be able to expand fast enough to maintain a population ahead of the advancing climate boundary. However, species that exhibit these traits are also more likely to become invasives themselves, which is the primary argument against assisted migration in the first place [36]. In this paper, we introduce a mathematical framework for studying species range expansion in the context of climate change. Because, however, this is a broad class of problems, with a number of different applications to relevant ecological questions, a variety of extensions are possible. First, it should be noted that the reactiondiffusion model we consider is a limiting case scenario. In particular, it allows for establishment of small populations at large distances that may not reflect the true dispersal behavior of species with finite maximum dispersal distances. While this is not expected to alter conclusions for species with c > c∗ (∞), species with c < c∗ (∞) may still go extinct, particularly if their maximum dispersal distances are small and range boundaries move too quickly. Extending our analysis to alternate model formulations, for example, integro-difference equations, would make it easier to capture this aspect of real systems. Integro-difference equations would also allow us to incorporate discrete annual life-histories that are common in certain classes of species, for example, insects in seasonal environments. Another interesting extension would be to consider the effects of alternate growth rate functions on the spread speed thresholds necessary to prevent extinction. In the current model, for example, we assume r = r(x − ct). However, r could also be periodic in space (e.g., reflecting habitat heterogeneity, as in [39]) or periodic in time (e.g., reflecting seasonality) and this could alter predictions regarding whether a species will go extinct. Spatially or temporally varying diffusion rates, d, could likewise be considered. The study of invasion dynamics in the context of climate change is in its infancy. For this paper, our goal was to outline a limiting case scenario that could be used to define the general behavior of a broad class of systems relevant to invasion dynamics under global change. In doing so, we showed that whether an expanding population can continue to spread in the face of ongoing change depends on how well the species disperses as well its intrinsic maximum rate of population growth. This helps us to identify which non-equilibrium populations are likely to continue their spread and which are likely to go extinct as a result of being overtaken by unsuitable habitat. In general, the species traits that lead to population persistence are those associated with species invasiveness. This leads us to conclude that climate change should have a minimal impact on current invasive species but a strong impact on native species still recovering from the most recent ice age. Our results also help to highlight the issues associated with assisted migration. In particular, introduced species that can spread at rates fast enough to outrun climate change are also those most likely to become invasives themselves.
SPREAD WITH A SHIFTING HABITAT EDGE
1415
Acknowledgments. The authors would like to thank two anonymous referees for their constructive comments which have considerably improved the paper. REFERENCES [1] D. A. Andow, P. M. Kareiva, S. A. Levin, and A. Okubo, Spread of invading organisms, Landscape Ecology, 4 (1990), pp. 177–188. [2] D. G. Aronson and H. Weinberger, Nonlinear diffusion in population genetics, combustion, and nerve impulse propagation, in Partial Differential Equations and Related Topics, Lecture Notes in Math. 446, Springer-Verlag, Berlin, 1975, pp. 5–49. [3] D. G. Aronson and H. F. Weinberger, Multidimensional nonlinear diffusion arising in population genetics, Adv. Math., 30 (1978), pp. 33-76. [4] M. B. Araujo and R. G. Pearson, Equilibrium of species’ distributions with climate, Ecogeography, 28 (2005), pp. 693–695. [5] H. Berestycki, O. Diekmann, C. J. Nagelkerke, and P. A. Zegeling, Can a species keep pace with a shifting climate?, Bull. Math. Biol., 71 (2009), pp. 399–429. [6] H. J. B. Birks and K. J. Willis, Alpines, trees, and refugia in Europe, Plant Ecology Diversity, 1 (2008), pp. 147–160. [7] B. A. Bradley, D. S. Wilcove, and M. Oppenheimer, Climate change increases risk of plant invasion in the eastern United States, Biological Invasions, 12 (2010), pp. 1899–1872. [8] J. R. Britton, J. Cucherousset, G. D. Davies, M. J. Godard, and G. H. Copp, Nonnative fishes and climate change: Predicting species responses to warming temperatures in temperate regions, Freshwater Biology, 55 (2010), pp. 1130–1141. [9] M. T. Burrows, D. S. Schoeman, L. B. Buckley, P. Moore, E. S. Poloczanska, K. M. Brander, C. Brown, J. F. Bruno, C. M. Duarte, B. S. Halpern, J. Holding, C. V. Kappe, W. Kiessling, M. I. O’Connor, J. M. Pandolfi, C. Paremsan, F. Schwing, W. J. Sydeman, and A. J. Richardson, The pace of shifting climate in marine and terrestria ecosystems, Science, 334 (2011), pp. 652–655. [10] W. F. Fagan, R. S. Cantrell, and C. Cosner, How habitat edges change species interactions, American Naturalist, 153 (1999), pp. 165–182. [11] W. F. Fagan, R. S. Cantrell, C. Cosner, and S. Ramakrishnan, Interspecific variation in critical patch size and gap crossing ability as determinants of geographic range size distributions, American Naturalist, 173 (2009), pp. 363–375. [12] R. A. Fisher, The wave of advance of advantageous genes, Ann. Eugenics, 7 (1937), pp. 355–369. [13] R. V. Gallagher, D. E. Duursma, J. O’Donnell, P. D. Wilson, P. O. Downey, L. Hughes, and M. R. Leishman, The grass may not always be greener: Projected reductions in climatic suitability for exotic grasses under future climates in Australia, Biological Invasions, 15 (2013), pp. 961–975. [14] P. Gonzalez, R. P. Neilson, J. M. Lenihan, and R. J. Drapek, Global patterns in the vulnerability of ecosystems to vegetation shifts due to climate change, Global Ecology Biogeography, 19 (2010), pp. 755–768. [15] J. J. Hellmann, J. E. Byers, B. G. Bierwagen, and J. S. Dukes, Five potential consequences of climate change for invasive species, Conservation Biology, 22 (2008), pp. 534–543. [16] J. Huang and W. Shen, Speeds of spread and propagation of KPP models in time almost and space periodic media, SIAM J. Appl. Dyn. Syst., 8 (2009), pp. 790–821. [17] J. M. Jeschke and D. L. Strayer, Usefulness of bioclimatic models for studying climate change and invasive species, Ann. New York Acad. Sci., 1134 (2008), pp. 1–24. ´ [18] A. Kolmogorov, I. Petrovskii, and N. Piscounov, Etude de l’´ equation de la diffusion avec croissance de la quantit´ e de mati´ ere et son application a un probl´ eme biologique, Moscou Univ. Math. Bull., 1 (1937), p. 126. [19] B. Li, H. F. Weinberger, and M. A. Lewis, Spreading speeds as slowest wave speeds for cooperative systems, Math. Biosci., 196 (2005), pp. 82–98. [20] S. R. Loarie, P. B. Duffy, H. Hamilton, G. P. Asner, C. B. Field, and D. D. Ackerly, The velocity of climate change, Nature, 462 (2009), pp. 1052–1055. [21] M. A. Lewis, B. Li, and H. F. Weinberger, Spreading speeds and the linear conjecture for two-species competition models, J. Math. Biol., 45 (2002), pp. 219–233. [22] H. J. Lynch, R. Naveen, P. N. Trathan, and W. F. Fagan, Environmental change and the shifting balance among penguins on the Antarctic Peninsula, Ecology, 93 (2012), pp. 1367–1377.
1416
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
[23] J. S. McLachlan, J. J. Hellmann, and M. W. Schwartz, A framework for debate of assisted migration in an era of climate change, Conservation Biology, 21 (2007), pp. 297–302. [24] L. W. Morrison, M. D. Korzukhin, and S. D. Porter, Predicted range expansion of the invasive fire ant, Solenopsis invicta, in the eastern United States based on the VEMAP global warming scenario, Diversity Distributions, 11 (2005), pp. 199–204. [25] J. Ni, A simulation of biomes on the Tibetan Plateau and their responses to global climate change, Mt. Res. Dev., 20 (2000), pp. 80–89. [26] S. Normand, R. E. Ricklefs, F. Skov, J. Bladt, O. Tackenberg, and J.-C. Svenning, Postglacial migration supplements climate in determining plant species ranges in Europe, Proc. Roy. Soc. B Biological Sci., 278 (2011), pp. 3644–3653. [27] C. V. Pao, Nonlinear Parabolic and Elliptic Equations, Plenum Press, New York, 1992. [28] A. Paradis, J. Elkinton, K. Hayhoe, and J. Buonaccorsi, Role of winter temperatures and climate change on the survival and future range expansion of the hemlock wooly adelgid (Adelges tsugae) in eastern North America, Mitigation Adaptation Strategies Global Change, 13 (2008), pp. 541–554. [29] F. Parker-Allie, C. F. Musil, and W. Thuiller, Effects of climate warming on the distributions of invasive Eurasion annual grasses: A South African perspective, Climate Change, 94 (2009), pp. 87–103. [30] C. Parmesan, Ecological and evolutionary responses to recent climate change, Ann. Rev. Ecology Evolution Systematics, 37 (2006), pp. 637–669. [31] C. L. Parr, E. F. Gray, and W. J. Bond, Cascading biodiversity and functional consequences of a global change-induced biome switch, Diversity Distributions, 18 (2012), pp. 493–503. [32] J. J. Polovina, J. P. Dunne, P. A. Woodworth, and E. A. Howell, Projected expansion of the subtropical biome and contraction of the temperate and equatorial upwelling biomes in the North Pacific under global warming, ICES J. Marine Sci., 68 (2011), pp. 986–995. [33] A. B. Potapov and M. A. Lewis, Climate and competition: The effect of moving range boundaries on habitat invisibility, Bull. Math. Biol., 66 (2004), pp. 975–1008. [34] F. J. Rahel and J. D. Olden, Asessing the effects of climate change on aquatic invasive species, Conservation Biology, 22 (2008), pp. 521–533. [35] M. Rejmanek, and D. M. Richardson, What Attributes Make Some Plant Species More Invasive? Ecology, 77 (1996), pp. 1655–1661. [36] A. Ricciardi and D. Simberloff, Assisted colonization is not a viable conservation strategy, Trends Ecology Evolution, 24 (2009), pp. 248–253. [37] M. Scheffer, M. Hirota, M. Holmgren, E. H. Van Nes, and F. S. Chapin, Thresholds for boreal biome transitions, Proc. Natl. Acad. Sci. USA, 109 (2012), pp. 21384–21389. [38] S. Scheiter and S. I. Higgins, Impacts of climate change on the vegetation of Africa: An adaptive dynamic vegetation modelling approach, Global Change Biology, 15 (2009), pp. 2224– 2246. [39] N. Shigesada, K. Kawasaki, and E. Teramoto, Traveling periodic waves in heterogeneous environments, Theoretical Population Biology, 30 (1986), pp. 143–160. [40] D. A. Smale and T. Wernberg, Extreme climatic event drives range contraction of a habitatforming species, Proc. Roy. Soc. B, 280 (2013). [41] R. W. Sutherst, Climate change and invasive species: A conceptual framework, in Invasive Species in a Changing World, H. A. Mooney, and R. J. Hobbs, eds., Island Press, Washington, DC, 2000, pp. 211–240. [42] J. C. Svenning and F. Skov, Ice age legacies in the geographical distribution of tree species richness in Europe, Global Ecology Biogeography, 16 (2007), pp. 234–245. [43] J. C. Svenning and F. Skov, Could the tree diversity pattern in Europe be generated by postglacial dispersal limitation? Ecology Letters, 10 (2007), pp. 453–460. [44] D. Y. P. Tng, B. P. Murphy, E. Weber, G. Sanders, G. J. Williamson, J. Kemp, and D. M. J. S. Bowman, Humid tropical rain forest has expanded into eucalypt forest and savanna over the last 50 years, Ecology Evolution, 2 (2012), pp. 34–45. [45] G. R. Walther, A. Roques, P. E. Hulme, M. T. Sykes, P. Pysek, I. Kuhn, M. Zobel, S. Bacher, Z. Botta-Dukat, H. Bugmann, B. Czucz, J. Dauber, T. Hickler, V. Jarosik, M. Kenis, S. Klotz, D. Minchin, M. Moora, W. Nentwig, J. Ott, V. Panov, B. Reineking, C. Robinet, V. Semenchenko, W. Solarz, W. Thuiller, M. Vila, K. Vohland, and J. Settele, Alien species in a warmer world: Risks and opportunities, Trends Ecology Evolution, 24 (2009), pp. 686–693. [46] X. Wang, On the Cauchy problem for reaction-diffusion equations, Trans. Amer. Math. Soc., 337 (1993), pp. 549–589.
SPREAD WITH A SHIFTING HABITAT EDGE
1417
[47] H. F. Weinberger, Long-time behavior of a class of biological models, SIAM J. Math. Anal., 13 (1982), pp. 353–396. [48] H. F. Weinberger, On spreading speeds and traveling waves for growth and migration models in a periodic habitat: Long-time behavior of a class of biological models, J. Math. Biol., 45 (2002), pp. 511–548. [49] H. F. Weinberger, M. A. Lewis, and B. Li, Analysis of linear determinacy for spread in cooperative models, J. Math. Biol., 45 (2002), pp. 183–218. [50] K. A. With, The landscape ecology of invasive spread, Conservation Biology, 16 (2002), pp. 1192–1203. [51] Y. Zhou and M. Kot, Discrete-time growth-dispersal models with shifting species ranges, Theoretical Ecology, 4 (2011), pp. 13–25.
c 2014 Society for Industrial and Applied Mathematics
PERSISTENCE AND SPREAD OF A SPECIES WITH A SHIFTING HABITAT EDGE∗ BINGTUAN LI† , SHARON BEWICK‡ , JIN SHANG† , AND WILLIAM F. FAGAN‡ Abstract. We study a reaction-diffusion model that describes the growth and spread of a species along a shifting habitat gradient on which the species’ growth increases. It is assumed that the linearized species growth rate is positive near positive infinity and is negative near negative infinity. We show that the persistence and spreading dynamics depend on the speed of the shifting habitat edge c and a number c∗ (∞) that is determined by the maximum linearized growth rate and the diffusion coefficient. We demonstrate that if c > c∗ (∞), then the species will become extinct in the habitat, and that if c < c∗ (∞), then the species will persist and spread along the shifting habitat gradient at an asymptotic spreading speed c∗ (∞). Key words. reaction-diffusion equation, shifting habitat edge, persistence, spreading speed AMS subject classifications. 92D40, 92D25 DOI. 10.1137/130938463
1. Introduction. Whether a species can persist and spread in an environment is a fundamental and long-standing question in ecology. Recently, this question has been brought to the forefront by anthropogenic disturbance, including climate warming and landscape conversion. In particular, anthropogenic change has been blamed for widespread population declines across a broad range of taxa [30]. However, at the same time, there is evidence that other populations, most notably pest species and alien invasives, may be benefitting from recent changes [7, 24, 28]. Why certain species suffer while others tolerate and even prosper under changing conditions remains an area of active investigation. Early on, it was recognized that species persistence and spread often depends on spatial context, specifically the combined effects of habitat range and species dispersal. This led to a vast body of literature analyzing the role of spatial processes in population dynamics [12, 18]. Many of the conclusions from these early studies pertain to species persistence in changing landscapes as well. Indeed, because drivers like climate change alter the spatial spread of suitable habitats, population level effects can often be understood in terms of perturbations to spatial processes. While several frameworks exist for studying spatial population ecology, much analysis is based on reaction-diffusion models like the Fisher’s equation (1.1)
∂2u ∂u = d 2 + ru − u2 . ∂t ∂x
Here u is the population density of the species of interest, d is the diffusion coefficient, r is the population growth rate, and −u2 accounts for density-dependent death (e.g., ∗ Received by the editors September 25, 2013; accepted for publication (in revised form) May 29, 2014; published electronically September 16, 2014. http://www.siam.org/journals/siap/74-5/93846.html † Department of Mathematics, University of Louisville, Louisville, KY 40059 (bing.li@louisville. edu, [email protected]). The research of these authors was partially supported by the National Science Foundation under grant DMS-1225693. ‡ Department of Biology, University of Maryland, College Park, MD 20742 (sharon bewick@ hotmail.com, [email protected]). The research of these authors was partially supported by the National Science Foundation under grant DMS-1225917.
1397
1398
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
due to resource limitation). Frequently, models such as (1.1) are used to study the spreading speed of a population. Loosely speaking, this is the asymptotic rate at which a species, initially introduced in a bounded domain, expands its spatial range. Ecologically, an understanding of spreading speeds can provide insight into invasion processes, for example, how quickly an introduced species can move into a novel landscape or how rapidly an extirpated species can recover to its previous √ range. It is well known that for the Fisher’s equation, the spreading speed is c∗ = 2 dr, which is the slowest speed of a class of traveling wave solutions to (1.1) [12, 18]. A traveling wave solution describes the propagation of a species as a wave with a fixed shape and a fixed speed. A number of extensions of the basic Fisher’s equation exist. It has, for example, been studied in systems with general growth functions [2], in multiple dimensions [3], and in systems with interacting species [10, 19, 21, 48]. Modified versions have been applied to discrete time systems in order to capture seasonal life-histories. Finally, a number of different forms of spatial heterogeneity have been considered [1, 11, 39, 48, 50]. To incorporate spatial heterogeneity, one or both of the diffusion coefficient, d, or the population growth rate, r, are taken to be functions of space. Shigesada, Kawasaki, and Teramoto, for example, studied spreading speeds in systems where both the diffusion coefficient and the growth rate are periodic [39]. Weinberger considered a more general problem for both continuous time and discrete time models [49]. In both papers, periodicity was attributed to patchiness in the environment. Spatial heterogeneity has also been considered for models in time almost periodic and space periodic media [16]. Recently, there has been renewed interested in spatial ecology, largely driven by the threats associated with global change. In particular, due to the processes of global change, habitat regions suitable for population growth of many species have shifted geographically and are predicted to shift even more. For example, changes in the boundaries between forest and non-forest habitat types are receiving extensive study by ecologists and biogeographers in biomes worldwide (e.g., [14, 37, 38, 44]). In many cases, the boundary between suitable and unsuitable habitat shifts gradually over a period of years or decades (e.g., [31, 32]), but in other cases the range edge can shift quickly [37] or even in response to a single climatic event [40]. Changes in habitat boundaries may result in shrinkage of particular habitat types, such as when forest habitat expands upslope, encroaching on alpine communities [14, 25]. Similarly, expansion of scrub-woodland habitats into areas historically dominated by grasslands decreases the spatial extent of the grasslands [14, 31]. In other situations, suitable habitat is actually expanding [14, 25]. For example, along the Antarctic Peninsula, temperatures have been increasing rapidly for several decades, resulting in a decline in spring sea ice. This has increased accessibility of ice-free territory on shore, allowing the gentoo penguin (Pygoscelis papua) to spread southward beyond its historical limits [22]. Understanding and predicting the effects of changes like those described above has required the development of new theoretical models. In general, the emphasis of these models is on the role of changing habitat suitability. Thus while the reactiondiffusion framework in (1.1) is still valid, it requires population growth rates that are functions of time and location. One of the most common effects of global change is shifts or translations in habitat ranges. This is often modeled by assuming that the population growth rate can be described as r(x − ct), where r(x) is the spatially varying baseline, or historic rate of population growth. Initially suitable habitat is all locations with r(x) > 0. As global change progresses, however, the spatial locations of
SPREAD WITH A SHIFTING HABITAT EDGE
1399
the suitable habitats move in the positive direction at a rate of c units per timestep. Berestycki et al. [5] used a model like this to define the minimum size of a moving habitat necessary to sustain a non-zero species population (analogous to the critical patch size in models with stationary habitats). They then generalized their analysis to derive criteria for the persistence of a species in any region with a moving and spatially varying habitat. In a similar study, Potapov and Lewis [33] examined the implications of a moving habitat on competition between two species. More recently, Zhou and Kot [51] have extended the work in Berestycki et al. [5] to consider discrete time systems. Today, most applications of diffusion equations to global change scenarios have focused on determining criteria for the persistence of an established species. These are species that exist at an equilibrium distribution at the onset of climate change but then must track climate change in order to persist into the future. Fewer studies have considered the spread or invasion of an introduced species and how this might be impacted by shifts in habitat suitability. In [33], the authors consider the invasion problem from the perspective of two competing species and a moving boundary. Here, we take a broader and more basic perspective and consider the spread of a single species over a region with varying habitat suitability that is shifting in time. This scenario is relevant to a number of different problems in ecology, for example, predicting the spread of an invasive species under climate change [15, 41] or ascertaining the likelihood of success of assisted migration (i.e., deliberate translocation of species to areas with more suitable conditions) [23]. Indeed, understanding how climate change will affect non-equilibrial populations like those associated with invasion processes is a broad class of problems—one that has garnered some attention from biological communities [15, 34, 41] but has been largely ignored from a mathematical modeling perspective (but see Jeschke and Strayer [17]). To explore the issue of species spread in the context of climate change we consider the following formulation of the reaction-diffusion model in (1.1): (1.2)
∂2u ∂u = d 2 + ur(x − ct) − u2 . ∂t ∂x
Specifically, we assume an infinite domain −∞ < x < +∞ and a constant diffusion coefficient, d. For population growth, we take r(x − ct), where c > 0 and r(ξ) is continuous and nondecreasing and bounded with r(−∞) < 0 and r(∞) > 0. r(x − ct) thus divides the spatial domain into two parts: the region with good-quality habitat suitable for growth (i.e., r(x − ct) > 0), and the region with poor-quality habitat unsuitable for growth (i.e., r(x − ct) < 0). The edge of the habitat suitable for species growth is shifting at a speed c. We show that the persistence and spreading dynamics depend on c and c∗ (∞) = 2 dr(∞). We demonstrate that if c > c∗ (∞), then the species will become extinct in the habitat, and that if c < c∗ (∞), then the species will persist and spread along the shifting habitat gradient at an asymptotic spreading speed c∗ (∞). This paper is organized as follows. In the next section, the mathematical results regarding the spatial dynamics of (1.2) are presented. Section 3 presents a number of simulations that help to illustrate the results from section 2. Finally, in section 4 we discuss the relevance and implications of our results. 2. Mathematical results. In this section we provide mathematical results on the spatial dynamics for model (1.2). We begin with making the following hypothesis for r(x).
1400
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
Hypothesis 2.1. r(x) is continuous, nondecreasing and bounded, and piecewise continuously differentiable in x for −∞ < x < ∞, 0 < r(∞) < ∞, and −∞ < r(−∞) < 0. For 0 ≤ u1 , u2 ≤ r(∞), −∞ < x < ∞, and t ≥ 0, |u1 (r(x − ct) − u1 ) − u2 (r(x − ct) − u2 )| ≤ 3r(∞)|u1 − u2 |, so that u(r(x − ct) − u) is Lipschitz continuous in u. Clearly u ≡ 0 is a trivial solution (thus a lower solution) and u ≡ r(∞) is an upper solution of (1.2). The theory on the existence and uniqueness of solutions for reaction-diffusion systems has been well established (e.g., [27, Theorem 2.1], [46, Lemma 1.2]). It is known that the initial value problem of (1.2) with u(0, x) = u0 (x) where u0 (x) is continuous and 0 ≤ u0 (x) ≤ r(∞) has a unique classical solution u(t, x) with 0 ≤ u(t, x) ≤ r(∞). Let ρ > 3r(∞). Then u(ρ + r(x − ct) − u) is nondecreasing in u for 0 ≤ u ≤ r(∞). Equation (1.2) can be written as ∂2u ∂u + ρu = d 2 + u(ρ + r(x − ct) − u). ∂t ∂x
(2.1)
The solution of (2.1) with u(0, x) = u0 (x) satisfies the integral equation +∞ u(t, x) = k(t, x − y)u0 (y)dy −∞ t
+∞
+ 0
−∞
k(t − τ, x − y)u(τ, y)[ρ + r(y − cτ ) − u(τ, y)]dydτ,
where y2 1 e−ρs− 4ds . k(s, y) = √ 4πds
Consider the sequence u(n) (t, x) defined by (2.2) u
(n+1)
(t, x) =
+∞
−∞ t
k(t, x − y)u0 (y)dy
+∞
+ 0
−∞
k(t − τ, x − y)u(n) (τ, y)[ρ + r(y − cτ ) − u(n) (τ, y)]dydτ,
where u(0) (t, x) = 0 or u(0) (t, x) = r(∞). Lemma 7.22 in [27] shows that if u(0) (t, x) = 0, then u(n) (t, x) is nondecreasing in n, and if u(0) (t, x) = r(∞), then u(n) (t, x) is nonincreasing in n, 0 ≤ u(n) (t, x) ≤ r(∞). In both cases, u(t, x) = limn→∞ u(n) (t, x) is the solution of (2.1) with u(0, x) = u0 (x), and 0 ≤ u(t, x) ≤ r(∞). We first provide a useful lemma for the equation (2.3)
∂2u ∂u = d 2 + u(r(x) − u). ∂t ∂x
Lemma 2.1. Let u ¯(t, x) be the solution of (2.3) with u ¯(0, x) = r(∞). Then u ¯(t, x) ¯(t, −∞) = 0, and u ¯(t, ∞) ≡ is nonincreasing in t and nondecreasing in x, limt→∞ u r(∞) for t > 0.
SPREAD WITH A SHIFTING HABITAT EDGE
1401
Proof. It follows from (2.2) with u(0) (t, x) = r(∞) and c = 0 that the solution u ¯(t, x) of (2.3) with u ¯(0, x) = r(∞) has the property that u ¯(t, x) = limn→∞ u(n) (t, x), where (2.4) u(n+1) (t, x) = e−ρt r(∞) t +∞ + k(s, y)u(n) (t − s, x − y)[ρ + r(x − y) − u(n) (t − s, y)]dyds. 0
−∞
Here we have used the simple fact that +∞ k(t, y)dy = e−ρt . −∞
We first consider n = 1 in (2.4). Direct calculations show that +∞ ∂u(1) (t, x) k(t, y)r(∞)[ρ + r(x − y) − r(∞)]dy = −ρe−ρt r(∞) + ∂t −∞ +∞ k(t, y)r(∞)[r(x − y) − r(∞)]dy = −∞
≤0 and ∂u(1) (t, x) = ∂x Assume that
∂u(n) (t,x) ∂t
t 0
≤ 0 and
+∞
−∞
k(s, y)r(∞)r (x − y)dyds ≥ 0.
∂u(n) (t,x) ∂x
≥ 0 for some n > 0. Then
∂u(n+1) (t, x) = −ρe−ρt r(∞) ∂t +∞ k(t, y)u(n) (0, x − y)[ρ + r(x − y) − u(n) (0, x − y)]dy + −∞
t
+∞
+ 0
−∞
k(s, y)
∂u(n) (t, x) [ρ + r(x − y) − 2u(n) (t − s, y)]dyds. ∂t
Note that ρ + r(x − y) − 2u(n) (t − s, y) ≥ 0 because of ρ > 3r(∞) and that 0 ≤ u(n) (t − s, y) ≤ r(∞). Since (2.4) shows u(n) (0, x) = r(∞) for all n, we have that +∞ ∂u(n+1) (t, x) = k(t, y)r(∞)[r(x − y) − r(∞)]dy ∂t −∞ t +∞ ∂u(n) (t, x) [ρ + r(x − y) − 2u(n) (t − s, y)]dyds k(s, y) + ∂t 0 −∞ ≤ 0. On the other hand, (n) t +∞ ∂u (t − s, x − y) ∂u(n+1) (t, x) = (ρ + r(x − y) − 2u(n) (t − s, x − y)) k(s, y) ∂x ∂x 0 −∞ (n) + r (x − y)u (t − s, x − y) dyds ≥ 0.
1402
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN (n)
By induction, ∂u limn→∞ u(n) (t, x),
(n)
(t,x) ∂t ∂u ¯ (t,x) ∂t
(2.5) u ¯(t, x) = r(∞)e−ρt +
≤ 0 and ∂u ∂x(t,x) ≥ 0 for all n > 0. Since u ¯(t, x) = ∂u ¯ (t,x) ≤ 0 and ∂x ≥ 0. u ¯ satisfies
t 0
+∞
−∞
k(s, y)¯ u(t − s, x − y)[ρ + r(x − y) − u ¯(t − s, x − y)]dyds.
Taking the limit x → −∞ in (2.5) and applying the dominated convergence theorem we obtain t +∞ ¯(t − s, −∞)] k(s, y)dyds u ¯(t, −∞) = r(∞)e−ρt + u¯(t − s, −∞)[ρ + r(−∞) − u 0
= r(∞)e−ρt + = r(∞)e−ρt +
−∞
t
0
0
t
e−ρs u ¯(t − s, −∞)[ρ + r(−∞) − u ¯(t − s, −∞)]ds e−ρ(t−s) u ¯(s, −∞)[ρ + r(−∞) − u ¯(s, −∞)]ds.
This implies that u ¯(t, −∞) satisfies the differential equation ∂u ¯(t, −∞) =u ¯(t, −∞)[r(−∞) − u ¯(t, −∞)]. ∂t Since r(−∞) < 0, 0 is the only nonnegative equilibrium of the differential equation and it attracts all nonnegative solutions. We therefore have limt→∞ u ¯(t, −∞) = 0. Similarly one can show that u ¯(t, ∞) satisfies the differential equation ∂u ¯(t, ∞) =u ¯(t, ∞)[r(∞) − u ¯(t, ∞)], ∂t which has two equilibria 0 and r(∞). Since u ¯(0, ∞) = r(∞), we have that u ¯(t, ∞) ≡ r(∞). The proof is complete. For r(x) > 0, define c∗ (x) = 2 dr(x). It is easily seen that c∗ (x) = inf φ(x; μ), μ>0
where φ(x; μ) = The infimum occurs at μ∗ (x) =
r(x) d .
dμ2 + r(x) . μ
The function
ψ(μ) = 2dμ is useful. It is easily seen that φ(x; μ) > ψ(μ) for 0 < μ < μ∗ (x) and φ(x; μ∗ (x)) = ψ(μ∗ (x)).
SPREAD WITH A SHIFTING HABITAT EDGE
1403
Our first theorem shows that if c∗ (∞) < c, then the species will eventually become extinct in space under certain conditions. Theorem 2.1. Assume that Hypothesis 2.1 is satisfied. Let c > 0 and c∗ (∞) < c. If 0 ≤ u0 (x) ≤ r(∞) and u0 (x) ≡ 0 for all sufficiently large x, then for every ε > 0 there exits T > 0 such that for t ≥ T , the solution u(t, x) of (1.2) with u(0, x) = u0 (x) satisfies u(t, x) ≤ ε for all x. Proof. Let u ¯(t, x) be the solution of (2.3) with u¯(0, x) = r(∞). Since ¯(t, −∞) = 0 by Lemma 2.1, for every ε > 0, there exist T1 > 0 and M > 0 limt→∞ u such that for x ≤ −M , u¯(T1 , x) < ε. Since u ¯(t, x) is nonincreasing in t, (2.6)
u ¯(t, x) < ε for t ≥ T1 and x ≤ −M.
The function u ˜(t, x) = u ¯(t, x − ct) satisfies (2.7)
∂u ˜ ∂2u ˜ ∂u ˜ =d 2 +c +u ˜(r(x − ct) − u ˜). ∂t ∂x ∂x
Since u ¯(t, x) is nondecreasing in x,
∂u ˜ ∂x
≥ 0. Since c > c∗ (∞) > 0, (2.7) shows that
∂ 2 u˜ ∂u ˜ ≥ d 2 + u˜(r(x − ct) − u ˜), ∂t ∂x so that u˜(t, x) = u¯(t, x − ct) with u˜(0, x) = r(∞) is an upper solution of (1.2). Since u(0, x) = u0 (x) ≤ r(∞), u(t, x) ≤ u ¯(t, x − ct). It follows from this and (2.6) that (2.8)
u(t, x) < ε for t ≥ T1 and x ≤ −M + ct.
Choose 0 < δ < c − c∗ (∞). Let μδ > 0 be the smaller positive solution of ∗ φ(∞; μ) = c∗ + δ/2. It is easily seen that u ˆ(t, x) = Ae−μδ (x−(c (∞)+δ/2)t) , where A is a positive constant, is a solution of the linear equation ∂2u ∂u = d 2 + r(∞)u. ∂t ∂x Since r(∞)u ≥ u(r(x − ct) − u), u ˆ(t, x) is an upper solution of (1.2). Choose A sufficiently large such that u0 (x) ≤ uˆ(0, x) = Ae−μδ x . Then u(t, x) ≤ Ae−μδ (x−(c
∗
(∞)+δ/2)t)
.
This shows that for x ≥ (c∗ (∞) + δ)t, u(t, x) ≤ Ae−(μδ δ/2)t . It follows that for the above given ε > 0, there exists T2 > 0 such that (2.9)
u(t, x) < ε for t ≥ T2 and x ≥ (c∗ (∞) + δ)t.
On the other hand, since c > c∗ (∞) + δ, there exists T3 > 0 such that for t ≥ T3 , −M +ct > (c∗ (∞)+δ)t. This, (2.8), and (2.9) show that for t ≥ T := max{T1 , T2 , T3 } and for all x, u(t, x) < ε. The proof is complete. Definition 2.2. We call a function u a continuous weak lower solution of model (1.2) if u is continuous for t ≥ 0 and −∞ < x < ∞, and ∂2u ∂u ≤ d 2 + u(r(x − ct) − u) ∂t ∂x
1404
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
in the distributional sense, i.e., for any T > 0 and any η ∈ C 2,1 ((−∞, ∞) × [0, T ]) with η ≥ 0 and supp η(·, t) being a bounded interval for all t ∈ [0, T ],
∞
−∞
≤
1 u(t, x)η(t, x)dx|t=T t=0
T1
∞
−∞
0
[u(s, x)(dηxx + ηt )(s, x) + η(s, x)u(s, x)(r(x − cs) − u(s, x))]dxds
if T1 ∈ [0, T ]. This definition is a slightly modified version of Definition 1.1 in [46]. Weak lower solutions were used in [2, 46] in studying reaction-diffusion systems. For γ > 0 and μ > 0, define (2.10)
v(μ; x) =
e−μx sin γx
if 0 ≤ x ≤ π/γ,
0
elsewhere.
(A scalar multiple of this function is denoted by v(s) in Weinberger [47].) v(μ; x) is a continuous function in x and its second order derivative in x exists and is continuous when x = 0, π/γ. The maximum of v(μ; x) occurs at σ(μ) = (1/γ) tan−1 (γ/μ). σ(μ) is a strictly decreasing function of μ. We have the following useful lemma. ∗
Lemma 2.3. Assume that c∗ (∞) > c ≥ 0. For any satisfying 0 < < c (∞)−c , 3 let be a number such that c∗ () = c∗ (∞) − . Let 0 < μ1 < μ2 < μ∗ () with ψ(μ1 ) = c + and ψ(μ2 ) = c∗ (∞) − 2. Then for any μ ∈ [μ1 , μ2 ], there exist a > 0 and γ > 0 sufficiently small such that av(μ; x − − ψ(μ)t) with v given by (2.10) is a continuous weak lower solution of (1.2). Furthermore, if u(0, x) ≥ av(μ; x − ), then u(t, x) ≥ av(μ; x − − ψ(μ)t) for all t > 0. Proof. By using the definition of v and integration by parts, we find that
T1
+∞
(2.11)
v(μ; x − − zs)(dηxx + ηt )(s, x)dxds
−∞ T1 +zs+π/γ
0
= 0
+zs T1 +zs+π/γ
0
+zs +zs+π/γ
=
+ +zs T1
+γ
0
e
v(μ; x − − zs)(dηxx + ηt )(s, x)dxds [dvxx (μ; x − − zs) − vs (μ; x − − zs)]η(s, x)dxds
v(μ; x − − zs)η(s, x)|T0 1 dxds
−μ π γ
π η s, + zs + + η(s, + zs) ds. γ
Direct calculations show that for μ ∈ [μ1 , μ2 ], x = + ψ(μ)t, and x = + ψ(μ)t + π/γ, (2.12)
dvxx (μ; x − − ψ(μ)t) − vt (μ; x − − ψ(μ)t) = −d(μ2 + γ 2 )v(μ; x − − ψ(μ)t).
SPREAD WITH A SHIFTING HABITAT EDGE
1405
Let f (t, x, u) = u(r(x − ct) − u). It follows from (2.11) and (2.12) that for μ ∈ [μ1 , μ2 ] and for sufficiently small a and γ, u ˜(μ; t, x) = av(μ; x − − ψ(μ)t) satisfies 0
T1
+∞
[˜ u(μ; s, x)(dηxx + ηt )(s, x) + η(s, x)f (s, x, u˜(μ; s, x))]dxds
−∞ ∞
−
=a
1 u˜(μ; s, x)η(s, x)dx|t=T t=0
−∞ T1 +ψ(μ)s+π/γ
0
[r(x − cs) − dμ2 − dγ 2 − av(μ; x − − ψ(μ)s)]
+ψ(μ)s
× v(μ; x − − ψ(μ)s)η(s, x)dxds T1 π −μ π γ + aγ e η s, + zs + + η(s, + zs) ds γ 0 T1 +ψ(μ)s+π/γ [r() − dμ22 − dγ 2 − av(μ; x − − ψ(μ)s)] ≥a 0
+ψ(μ)s
× v(μ; x − − ψ(μ)s)η(s, x)dxds T1 π −μ π γ η s, + zs + + aγ e + η(s, + zs) ds γ 0 T1 +ψ(μ)s+π/γ
∗ (c (∞) − 1.5) − dγ 2 − av(μ; x − − ψ(μ)s) =a 2d +ψ(μ)s 0 × v(μ; x − − ψ(μ)s)η(s, x)dxds T1 π −μ π γ + aγ η s, + zs + e + η(s, + zs) ds γ 0 T1 +ψ(μ)s+π/γ (c + ) − dγ 2 − av(μ; x − − ψ(μ)s) ≥a 2d +ψ(μ)s 0 × v(μ; x − − ψ(μ)s)η(s, x)dxds T1 π −μ π γ η s, + zs + + aγ e + η(s, + zs) ds γ 0 ≥ 0. It follows from Definition 2.2 that for μ ∈ [μ1 , μ2 ] and for sufficiently small a > 0 and γ > 0, av(μ; x − − ψ(μ)t) is a continuous weak lower solution of (1.2). If u(0, x) ≥ av(μ; x − ), then it follows from Lemma 1.2 in Wang [46] that u(t, x) ≥ av(μ; x − − ψ(μ)t) for all t > 0. The proof is complete. The following theorem shows that if c∗ (∞) > c, then the species persists in space and spreads to the right at the asymptotic spreading speed c∗ (∞). Theorem 2.2. Assume that Hypothesis 2.1 is satisfied. Let c∗ (∞) > c ≥ 0. Then the following statements are valid: (i) If 0 ≤ u(0, x) ≤ r(∞), then for any ε > 0, lim sup u(t, x) = 0. t→∞
x≤t(c−ε)
(ii) If 0 ≤ u(0, x) ≤ r(∞), and u(0, x) ≡ 0 for all sufficiently large x, then for any ε > 0,
1406
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
lim
t→∞
sup x≥t(c∗ (∞)+ε)
u(t, x) = 0.
(iii) If 0 ≤ u(0, x) ≤ r(∞), and u(0, x) > 0 on a closed interval, then for every ε with 0 < ε < (c∗ (∞) − c)/2, |r(∞) − u(t, x)| = 0. sup lim t→∞
t(c+ε)≤x≤t(c∗ (∞)−ε)
Proof. The first part of the proof of Theorem 2.1 shows that u ¯(t, x − ct) with u ¯(0, x) = r(∞) is an upper solution of (1.2). Since u(0, x) ≤ r(∞), u(t, x) ≤ u ¯(t, x − ct). Lemma 2.1 shows that for any > 0, there exist T1 > 0 and a real number M such that for t ≥ T1 and x ≤ −M + ct, u(t, x) < . On the other hand, for any ε > 0 there exists T2 > 0 such that for t ≥ T2 , (c − ε)t ≤ −M + ct. We therefore have that for t ≥ max{T1 , T2 } and x ≤ t(c − ε), u(t, x) < . This proves statement (i). The proof of statement (ii) is similar to the second half of the proof of Theorem 2.1 and is omitted. We now prove statement (iii). For any small positive with 0 < < ∗ }, let be a number such that min{ 13 , r(∞), c (∞)−c 3 c∗ () = c∗ (∞) − . Choose 0 < μ1 < μ2 < μ∗ () such that ψ(μ1 ) = c + and ψ(μ2 ) = c∗ (∞) − 2. By α v(μ; x−−ψ(μ)t) Lemma 2.3, for any μ ∈ [μ1 , μ2 ] and small α > 0 and γ > 0, v(μ;σ(μ)) with v given by (2.10) is a continuous weak lower solution of (1.2). Since u(0, x) ≥ 0 and u(0, x) ≡ 0, u(t, x) > 0 for all x and any t > 0. Choose σ(μ1 ) 0 < t0 < ψ(μ , and choose α and γ sufficiently small such that u(t0 , x) ≥ α for 1) x ∈ [, + 4π/γ]. Define ⎧ α ⎪ v(μ1 ; x − ) if ≤ x ≤ + σ(μ1 ), ⎪ ⎪ ⎪ v(μ ; 1 σ(μ1 )) ⎪ ⎪ ⎪ ⎨ α if + σ(μ1 ) ≤ x ≤ + 3π/γ + σ(μ2 ), w(0, x) = α ⎪ ⎪ v(μ2 ; x − − 3π/γ) if + 3π/γ + σ(μ2 ) ≤ x ≤ + 4π/γ, ⎪ ⎪ v(μ2 ; σ(μ2 )) ⎪ ⎪ ⎪ ⎩ 0 elsewhere. α v(μ1 ; x − − s) for 0 ≤ s ≤ 2π/γ, and It is easily seen that w(0, x) ≥ v(μ1 ;σ(μ 1 )) α w(0, x) ≥ v(μ2 ;σ(μ2 )) v(μ2 ; x − − 3π/γ + s) for 0 ≤ s ≤ 2π/γ. Since u(t0 , x) ≥ α for x ∈ [, + 4π/γ], Lemma 2.2 shows that for t ≥ t0 and 0 ≤ s ≤ 2π/γ,
(2.13)
u(t, x) ≥
α v(μ1 ; x − − ψ(μ1 )(t − t0 ) − s) v(μ1 ; σ(μ1 ))
SPREAD WITH A SHIFTING HABITAT EDGE
1407
and (2.14)
u(t, x) ≥
α v(μ2 ; x − − 3π/γ − ψ(μ2 )(t − t0 ) + s). v(μ2 ; σ(μ2 ))
Equation (2.13) implies that for t ≥ t0 , (2.15) u(t, x) ⎧ α ⎪ if + ψ(μ1 )(t − t0 ) ≤ x ≤ + σ(μ1 ) ⎪ ⎪ v(μ1 ; σ(μ1 )) ⎪ ⎪ ⎪ ⎪ ⎪ + ψ(μ1 )(t − t0 ), ⎨ × v(μ1 ; x − − ψ(μ1 )(t − t0 )) ≥ α if + σ(μ1 ) + ψ(μ1 )(t − t0 ) ≤ x ≤ + 2π/γ ⎪ ⎪ ⎪ ⎪ ⎪ + σ(μ1 ) + ψ(μ1 )(t − t0 ), ⎪ ⎪ ⎪ ⎩ 0 elsewhere. On the other hand, (2.14) indicates that for t ≥ t0 , (2.16) u(t, x) ⎧ α ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ α ≥ v(μ ; ⎪ 2 σ(μ2 )) ⎪ ⎪ ⎪ ⎪ × v(μ2 ; x − − 3π/γ − ψ(μ2 )(t − t0 )) ⎪ ⎪ ⎪ ⎩ 0
if + π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) ≤ x ≤ + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ), if + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) ≤ x ≤ + 4π/γ + ψ(μ2 )(t − t0 ), elsewhere.
1 )−σ(μ2 ) Let h = π/γ+σ(μ ψ(μ2 )−ψ(μ1 ) . Since + 2π/γ + σ(μ1 ) + ψ(μ1 )(t − t0 ) ≥ + π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) for t0 ≤ t ≤ t0 + h, (2.15) and (2.16) show that for t0 ≤ t ≤ t0 + h,
(2.17)
u(t, x) ≥ w(t − t0 , x),
where (2.18) w(t − t0 , x) ⎧ α ⎪ ⎪ ⎪ ⎪ v(μ ; 1 σ(μ1 )) ⎪ ⎪ ⎪ ⎪ ⎪ × v(μ1 ; x − − ψ(μ1 )(t − t0 )) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ α = ⎪ ⎪ α ⎪ ⎪ ⎪ ⎪ ⎪ v(μ ; 2 σ(μ2 )) ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ × v(μ2 ; x − − 3π/γ − ψ(μ2 )(t − t0 )) ⎪ ⎪ ⎩ 0
if + ψ(μ1 )(t − t0 ) ≤ x ≤ + σ(μ1 ) + ψ(μ1 )(t − t0 ), if + σ(μ1 ) + ψ(μ1 )(t − t0 ) ≤ x ≤ + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ), if + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) ≤ x ≤ + 4π/γ + ψ(μ2 )(t − t0 ), elsewhere.
1408
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
We claim that (2.17) is valid for all t ≥ t0 . Assume that (2.17) is true for t0 ≤ α t ≤ t0 + nh for some positive integer n. Then w(nh, x) ≥ v(μ1 ;σ(μ v(μ1 ; x − − 1 )) α v(μ ; x − − 3π/γ − nhψ(μ ) nhψ(μ1 ) − s) and w(nh, x) ≥ v(μ2 ;σ(μ 2 2 + s) for 0 ≤ 2 )) s ≤ 2π/γ + (ψ(μ2 ) − ψ(μ1 ))nh. It follows from Lemma 2.3 that for t ≥ t0 + nh and 0 ≤ s ≤ 2π/γ + (ψ(μ2 ) − ψ(μ1 ))nh, (2.19)
u(t, x) ≥
α v(μ1 ; x − − nhψ(μ1 ) − ψ(μ1 )(t − (t0 + nh)) − s) v(μ1 ; σ(μ1 ))
and (2.20) u(t, x) ≥
α v(μ2 ; x − − 3π/γ − nhψ(μ2 ) − ψ(μ2 )(t − (t0 + nh)) + s). v(μ2 ; σ(μ2 ))
Equation (2.19) shows that u(t, x) ≥ α
(2.21) for
t ≥ t0 + nh
(2.22) and (2.23)
+ σ(μ1 ) + nhψ(μ1 ) + ψ(μ1 )(t − (t0 + nh)) ≤ x ≤ + 2π/γ + σ(μ1 ) + nhψ(μ2 ) + ψ(μ1 )(t − (t0 + nh)).
Equation (2.20) implies that (2.21) holds if (2.22) is satisfied and if (2.24)
+ π/γ + σ(μ2 ) + nhψ(μ1 ) + ψ(μ2 )(t − (t0 + nh)) ≤ x ≤ + 3π/γ + σ(μ2 ) + nhψ(μ2 ) + ψ(μ2 )(t − (t0 + nh)).
The intervals described by (2.23) and (2.24) overlap if + 2π/γ + σ(μ1 ) + nhψ(μ2 ) + ψ(μ1 )(t − (t0 + nh)) ≥ + π/γ + σ(μ2 ) + nhψ(μ1 ) + ψ(μ2 )(t − (t0 + nh)). This is equivalent to t ≤ t0 + 2nh +
π/γ + σ(μ1 ) − σ(μ2 ) = t0 + (2n + 1)h. ψ(μ2 ) − ψ(μ1 )
It follows that the intervals described by (2.23) and (2.24) overlap if t0 + nh ≤ t ≤ t0 + (n + 1)h. We therefore have that (2.21) holds if t0 + nh ≤ t ≤ t0 + (n + 1)h and + σ(μ1 ) + nhψ(μ1 ) + ψ(μ1 )(t − (t0 + nh)) ≤ x ≤ + 3π/γ + σ(μ2 ) + nhψ(μ2 ) + ψ(μ2 )(t − (t0 + nh)). It follows from this, (2.19), and (2.20) that (2.17) holds for t0 + nh ≤ t ≤ t0 + (n + 1)h. By induction, (2.17) holds for all t ≥ t0 . For the chosen > 0, there is L > 0 such that L 2 1 √ e−x dx ≥ 1 − . π −L
SPREAD WITH A SHIFTING HABITAT EDGE
1409
Note that for any s > 0, √ L 4ds
√ −L 4ds
x2 1 √ e− 4ds dx = 4πds
L
−L
2 1 √ e−ξ dξ. π
Let t1 > t0 be a sufficiently large number. For t > t1 , the solution u(t, x) satisfies the integral equation +∞ (2.25) k(t − t1 , x − y)u(t1 , y)dy u(t, x) = −∞ t
+∞
+ t1
−∞
k(t − τ, x − y)u(τ, y)[ρ + r(y − cτ ) − u(τ, y)]dydτ.
It follows from this and (2.17) that for t > t1 , (2.26)
u(t, x) ≥
+∞
−∞ t
k(t − t1 , x − y)w(t1 − t0 , y)dy
+∞
+ −∞
t1
k(t − τ, x − y)w(τ − t0 , y)[ρ + r(y − cτ ) − w(τ − t0 , y)]dydτ.
For t > t1 , x satisfying (2.27)
+ σ(μ1 ) + ψ(μ1 )(t − t0 ) + L 4d(t − t1 ) ≤ x ≤ + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) − L
and y satisfying
4d(t − t1 ),
−L 4d(t − t1 ) ≤ y ≤ L 4d(t − t1 ),
(2.28) we have that
+ σ(μ1 ) + ψ(μ1 )(t − t0 ) ≤ x − y ≤ + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) and x − y − ct ≥ + σ(μ1 ) + ψ(μ1 )(t − t0 ) − ct = + t + σ(μ1 ) − ψ(μ1 )t0 > . It follows from this and (2.18) that for x satisfying (2.27), +∞ (2.29) k(t − t1 , x − y)w(t1 − t0 , y)dy −∞
+∞
k(t − t1 , y)w(t1 − t0 , x − y)dy L√4d(t−t1 ) y2 1 − −ρ(t−t1 ) ≥e e 4d(t−t1 ) w(t1 − t0 , x − y)dy √ 4πd(t − t1 ) −L 4d(t−t1 ) √ L 4d(t−t1 ) y2 1 − e 4d(t−t1 ) dy ≥ αe−ρ(t−t1 ) √ 4πd(t − t1 ) −L 4d(t−t1 ) =
−∞
≥ (1 − )αe−ρ(t−t1 )
1410
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
and
+∞
(2.30) −∞
k(t − τ, x − y)w(τ − t0 , y)[ρ + r(y − cτ ) − w(τ − t0 , y)]dy
≥ (1 − )α[ρ + r(∞) − − α]e−ρ(t−τ ) . Here we have used the fact that for x satisfying (2.27) and y satisfying (2.28) r(x − y − ct) ≥ r() > r(∞) − . It follows from (2.26), (2.29), and (2.30) that for t ≥ t1 and x satisfying (2.27) u(t, x) ≥ u ˜(1) (t), where u ˜
(1)
(t) = (1 − )αe
−ρ(t−t1 )
+ (1 − )
t
e−ρ(t−τ ) α[ρ + r(∞) − − α]dτ.
t1
Equation (2.25) and induction show that for t ≥ t1 and x satisfying + σ(μ1 ) + ψ(μ1 )(t − t0 ) + nL 4d(t − t1 ) (2.31) ≤ x ≤ + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) − nL 4d(t − t1 ), u(t, x) ≥ u ˜(n) (t),
(2.32) where u ˜(n) (t) satisfies (2.33)
u ˜(n) (t) = (1 − )αe−ρ(t−t1 ) t + (1 − ) e−ρ(t−τ ) u ˜(n−1) (τ )[ρ + r(∞) − − u˜(n−1) (τ )]dτ. t1
Direct calculations and induction show that (2.34)
u˜(n) (t) = an + bn (t)e−ρ(t−t1 ) ,
where an = (1 − )an−1 (ρ + r(∞) − − an−1 )/ρ, a1 = (1 − )α[ρ + r(∞) − − α]/ρ, and bn (t) is a sum of polynomials, and products of polynomials and exponential functions ˜(n) (t) = an . in the form of e−jρ(t−t1 ) with j a positive integer. Observe that limt→∞ u For small and α, an increases to r(∞) − − ρ/(1 − ) as n → ∞. This and (2.34) show that there exist a positive integer N and t2 > t1 such that for t > t2 and x satisfies (2.31) with n replaced by N , (2.35)
u ˜(N ) (t) ≥ r(∞) − 2 − ρ/(1 − ).
We choose t1 sufficiently large such that for t ≥ t1 , + σ(μ1 ) + ψ(μ1 )(t − t0 ) + N L 4d(t − t1 )
< + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) − N L 4d(t − t1 ).
For any given ε with 0 < ε < (c∗ (∞) − c)/2, choose sufficiently small such that < ε/2. Then there exists t3 > t2 such that for t > t3 , + σ(μ1 ) + ψ(μ1 )(t − t0 ) + N L 4d(t − t1 ) < t(c + ε) < t(c∗ (∞) − ε)
1411
SPREAD WITH A SHIFTING HABITAT EDGE
and t(c∗ (∞) − ε) < + 3π/γ + σ(μ2 ) + ψ(μ2 )(t − t0 ) − N L 4d(t − t1 ). This shows that for t > t3 , t(c + ε) ≤ x ≤ t(c∗ (∞) − ε) implies that x satisfies (2.31) with n replaced by N . It follows from (2.32) and (2.35) that u(t, x) ≥ r(∞) − 2 − ρ/(1 − ). (2.36) lim inf ∗ t→∞ t(c+ε)≤x≤t(c (∞)−ε)
Since is arbitrarily small and u(t, x) ≤ r(∞) for all x and t, (2.36) shows that statement (iii) holds. The proof is complete. 3. Simulations. In this section we present some numerical simulations to the model (1.2) with −0.5 if x ≤ ct, r(x − ct) = 1 elsewhere, where c > 0, and the initial data 0.5 sin(x − 20) if 20 ≤ x ≤ 20 + π, (3.1) u(0, x) = 0 elsewhere. r(x − ct) moves to the right at a speed c, and the species is initially introduced on the interval [20, 20 + π]. Choose d = 1 so that c∗ (∞) = 2. Numerical simulations were conducted using MATLAB. Figure 1 displays the numerical solution with c = 2.1 > c∗ (∞), which shows that the solution eventually becomes extinct. Figure 2 displays
1.2 t=35
Population density
1
t=90 0.8
0.6
t=1
0.4
0.2 t=165 0
0
50
100
150
200 250 Location(x)
300
350
400
Fig. 1. A numerical approximation to the graph of u(t, x) with c = 2.1 and the initial data given by (3.1). The solution finally becomes extinct.
1412
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
1.2 t=35
Population density
1
t=90
t=165
0.8
0.6 t=1 0.4
0.2
0
0
50
100
150
200 250 Location(x)
300
350
400
Fig. 2. A numerical approximation to the graph of u(t, x) with c = 1.5 and the initial data given by (3.1). The solution persists and spreads to the right at the asymptotic speed c∗ (∞) = 2.
the numerical solution with c = 1.5 < c∗ (∞), which shows that the solution expends its spatial range to the right at the asymptotic speed c∗ (∞) = 2. The numerical results support Theorems 2.1 and 2.2. 4. Discussion. The ability to persist in a changing landscape is one of the biggest challenges faced by species today. Indeed, because of global change, species are being forced to shift their historical ranges in order to keep up with shifts in suitable habitat. Most obviously, this involves movement poleward or to higher elevations. However, it can also include movement along moisture or precipitation gradients, salinity gradients, or other abiotic gradients impacted by global change. While the issue of persistence in a changing landscape has been addressed by a number of theoretical studies, these have primarily focused on established species that exist at equilibrium distributions in a bounded domain prior to the onset of climate change. In this case, species are restricted to and distributed over a band of suitable habitat. Persistence then requires that a species keep pace with the movement of its habitat band. This is a good model for species that are limited by abiotic constraints at both an upper and a lower boundary, as would be the case for native species restricted to a range of elevations or latitudes. A different scenario emerges for species that are established over a finite range but, at least for the forseeable future, only limited at one boundary. For this scenario, persistence requires that the species spread away from the boundary fast enough to outrun boundary encroachment. In contrast to the persistence of established species, this problem has received less attention. It is, however, a good model for invasive species that are initially introduced to a small portion of their suitable habitat. It is also a good model for native species that, because of some historical event, have not reached an equilibrium distribution across their native range. Examples include
SPREAD WITH A SHIFTING HABITAT EDGE
1413
certain tree species [42, 43], as well as a number of amphibian and reptile species [4] whose current distributions still show evidence of the Last Glacial Maximum (LGM). Finally, understanding when a species can outrun a moving boundary has implications for assisted migration [23]. If, for example, a species is moved from a shrinking habitat at the top of a mountain to a much larger habitat at a higher latitude, it would be useful to know whether the introduced population will persist, given the encroachment of unsuitable conditions from the south. Very generally, we find that a spreading species population will go extinct if its habitat boundary moves at speeds greater than the maximum rate of expansion of the species population, c > c∗ (∞). On the other hand, the population will persist if its habitat boundary moves at lower speeds c < c∗ (∞). Ultimately, this defines whether a species will be able to outrun an encroaching boundary. First, let us consider our results in the context of invasive species. A large number of invasives are currently cold-limited and thus are expected to benefit from climate change [7, 24, 28]. However, others will be negatively impacted. These species will lose habitat as a result of climate change, which may, in turn, limit their ability to persist and spread as alien invasives. This has led some authors to speculate that climate change will be a zero-sum game for invasive species [15]. In particular, while climate change will likely lead to the worsening or new emergence of a handful of invasives, it should reduce the impact of others, ultimately resulting in no change to the net burden of invasive species. Not surprisingly, less is known about species that show range contraction or reduced impacts as a result of climate change [45]. Nevertheless, several examples exist, for instance, the two fish species, sunbleak and ide, in England [8] as well as a number of Eurasion annual grasses in Southern Africa [29] and perennial grasses in Australia [13]. A key finding of our model is that despite losing habitat at one boundary, introduced species may continue spreading, provided that they can expand at a rate faster than climate change. In our model, the rate of expansion of a particular population depends on how well the species disperses, d, as well as its linearized maximum rate of population growth, r(∞). Interestingly, the majority of invasive species are characterized by rapid population growth and enhanced dispersal. In tree species, for instance, invasiveness has been associated with short juvenile periods, short intervals between large seed crops, and small seed size [35]. The first two traits lead to increased rates of population growth, while small seed size can increase dispersal. Consequently, most invasive species will be characterized by both large d and large r(∞), suggesting that climate change will rarely, if ever, halt the spread of an invasive species. This appears to be the case for some of the annual grasses in Southern Africa that are continuing their spread to high alpine communities with devastating effects [29]. The situation is different for native species with nonequilibrium distributions. This can be seen most dramatically with the example of species still spreading northward from ice-age refugia. Quite obviously, these species are expanding at a minimal rate. If this is because they are intrinsically poor at breeding and dispersal, then given estimates of current climate velocity on the order of 80 to 2700 m/yr [9, 20], it seems likely that species still strongly associated with ice-age refugia will have c > c∗ (∞). Consequently, existing populations of these species will not expand fast enough to escape encroachment by warmer climates. For northern hemisphere species with a southern range limit, the end result will be extinction. While a number of European species limited to ice-age refugia are Mediterranean and thus less likely to be severely limited by warmer climates, others are associated with alpine or Pannonian climate zones [26]. Species in these latter categories, for example, Pritzelago alpina or Pulsatilla alpina, are thought to have undergone range shifts and range
1414
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
contractions at their southern and/or low elevation boundaries in response to postglacial warming following the LGM [6, 26]. Consequently, these species are also likely to exhibit range losses in response to current warming regimes. Because they are limited by climate at their southern (low elevation) range boundary and by spreading speeds at their northern (high elevation) range boundary, we predict that alpine and Pannionian species that currently exhibit a strong association with ice-age refugia will be at a particularly high risk for extinction under ongoing global change. Moreover, transporting these species to higher latitudes will be of little use since, even there, they will ultimately be overtaken by climates that are unsuitable for growth. While the assisted migration of species from ice-age refugia may be futile, our model does not imply that all assisted migration will fail. Species with high dispersal and reproductive rates should be able to expand fast enough to maintain a population ahead of the advancing climate boundary. However, species that exhibit these traits are also more likely to become invasives themselves, which is the primary argument against assisted migration in the first place [36]. In this paper, we introduce a mathematical framework for studying species range expansion in the context of climate change. Because, however, this is a broad class of problems, with a number of different applications to relevant ecological questions, a variety of extensions are possible. First, it should be noted that the reactiondiffusion model we consider is a limiting case scenario. In particular, it allows for establishment of small populations at large distances that may not reflect the true dispersal behavior of species with finite maximum dispersal distances. While this is not expected to alter conclusions for species with c > c∗ (∞), species with c < c∗ (∞) may still go extinct, particularly if their maximum dispersal distances are small and range boundaries move too quickly. Extending our analysis to alternate model formulations, for example, integro-difference equations, would make it easier to capture this aspect of real systems. Integro-difference equations would also allow us to incorporate discrete annual life-histories that are common in certain classes of species, for example, insects in seasonal environments. Another interesting extension would be to consider the effects of alternate growth rate functions on the spread speed thresholds necessary to prevent extinction. In the current model, for example, we assume r = r(x − ct). However, r could also be periodic in space (e.g., reflecting habitat heterogeneity, as in [39]) or periodic in time (e.g., reflecting seasonality) and this could alter predictions regarding whether a species will go extinct. Spatially or temporally varying diffusion rates, d, could likewise be considered. The study of invasion dynamics in the context of climate change is in its infancy. For this paper, our goal was to outline a limiting case scenario that could be used to define the general behavior of a broad class of systems relevant to invasion dynamics under global change. In doing so, we showed that whether an expanding population can continue to spread in the face of ongoing change depends on how well the species disperses as well its intrinsic maximum rate of population growth. This helps us to identify which non-equilibrium populations are likely to continue their spread and which are likely to go extinct as a result of being overtaken by unsuitable habitat. In general, the species traits that lead to population persistence are those associated with species invasiveness. This leads us to conclude that climate change should have a minimal impact on current invasive species but a strong impact on native species still recovering from the most recent ice age. Our results also help to highlight the issues associated with assisted migration. In particular, introduced species that can spread at rates fast enough to outrun climate change are also those most likely to become invasives themselves.
SPREAD WITH A SHIFTING HABITAT EDGE
1415
Acknowledgments. The authors would like to thank two anonymous referees for their constructive comments which have considerably improved the paper. REFERENCES [1] D. A. Andow, P. M. Kareiva, S. A. Levin, and A. Okubo, Spread of invading organisms, Landscape Ecology, 4 (1990), pp. 177–188. [2] D. G. Aronson and H. Weinberger, Nonlinear diffusion in population genetics, combustion, and nerve impulse propagation, in Partial Differential Equations and Related Topics, Lecture Notes in Math. 446, Springer-Verlag, Berlin, 1975, pp. 5–49. [3] D. G. Aronson and H. F. Weinberger, Multidimensional nonlinear diffusion arising in population genetics, Adv. Math., 30 (1978), pp. 33-76. [4] M. B. Araujo and R. G. Pearson, Equilibrium of species’ distributions with climate, Ecogeography, 28 (2005), pp. 693–695. [5] H. Berestycki, O. Diekmann, C. J. Nagelkerke, and P. A. Zegeling, Can a species keep pace with a shifting climate?, Bull. Math. Biol., 71 (2009), pp. 399–429. [6] H. J. B. Birks and K. J. Willis, Alpines, trees, and refugia in Europe, Plant Ecology Diversity, 1 (2008), pp. 147–160. [7] B. A. Bradley, D. S. Wilcove, and M. Oppenheimer, Climate change increases risk of plant invasion in the eastern United States, Biological Invasions, 12 (2010), pp. 1899–1872. [8] J. R. Britton, J. Cucherousset, G. D. Davies, M. J. Godard, and G. H. Copp, Nonnative fishes and climate change: Predicting species responses to warming temperatures in temperate regions, Freshwater Biology, 55 (2010), pp. 1130–1141. [9] M. T. Burrows, D. S. Schoeman, L. B. Buckley, P. Moore, E. S. Poloczanska, K. M. Brander, C. Brown, J. F. Bruno, C. M. Duarte, B. S. Halpern, J. Holding, C. V. Kappe, W. Kiessling, M. I. O’Connor, J. M. Pandolfi, C. Paremsan, F. Schwing, W. J. Sydeman, and A. J. Richardson, The pace of shifting climate in marine and terrestria ecosystems, Science, 334 (2011), pp. 652–655. [10] W. F. Fagan, R. S. Cantrell, and C. Cosner, How habitat edges change species interactions, American Naturalist, 153 (1999), pp. 165–182. [11] W. F. Fagan, R. S. Cantrell, C. Cosner, and S. Ramakrishnan, Interspecific variation in critical patch size and gap crossing ability as determinants of geographic range size distributions, American Naturalist, 173 (2009), pp. 363–375. [12] R. A. Fisher, The wave of advance of advantageous genes, Ann. Eugenics, 7 (1937), pp. 355–369. [13] R. V. Gallagher, D. E. Duursma, J. O’Donnell, P. D. Wilson, P. O. Downey, L. Hughes, and M. R. Leishman, The grass may not always be greener: Projected reductions in climatic suitability for exotic grasses under future climates in Australia, Biological Invasions, 15 (2013), pp. 961–975. [14] P. Gonzalez, R. P. Neilson, J. M. Lenihan, and R. J. Drapek, Global patterns in the vulnerability of ecosystems to vegetation shifts due to climate change, Global Ecology Biogeography, 19 (2010), pp. 755–768. [15] J. J. Hellmann, J. E. Byers, B. G. Bierwagen, and J. S. Dukes, Five potential consequences of climate change for invasive species, Conservation Biology, 22 (2008), pp. 534–543. [16] J. Huang and W. Shen, Speeds of spread and propagation of KPP models in time almost and space periodic media, SIAM J. Appl. Dyn. Syst., 8 (2009), pp. 790–821. [17] J. M. Jeschke and D. L. Strayer, Usefulness of bioclimatic models for studying climate change and invasive species, Ann. New York Acad. Sci., 1134 (2008), pp. 1–24. ´ [18] A. Kolmogorov, I. Petrovskii, and N. Piscounov, Etude de l’´ equation de la diffusion avec croissance de la quantit´ e de mati´ ere et son application a un probl´ eme biologique, Moscou Univ. Math. Bull., 1 (1937), p. 126. [19] B. Li, H. F. Weinberger, and M. A. Lewis, Spreading speeds as slowest wave speeds for cooperative systems, Math. Biosci., 196 (2005), pp. 82–98. [20] S. R. Loarie, P. B. Duffy, H. Hamilton, G. P. Asner, C. B. Field, and D. D. Ackerly, The velocity of climate change, Nature, 462 (2009), pp. 1052–1055. [21] M. A. Lewis, B. Li, and H. F. Weinberger, Spreading speeds and the linear conjecture for two-species competition models, J. Math. Biol., 45 (2002), pp. 219–233. [22] H. J. Lynch, R. Naveen, P. N. Trathan, and W. F. Fagan, Environmental change and the shifting balance among penguins on the Antarctic Peninsula, Ecology, 93 (2012), pp. 1367–1377.
1416
B. LI, S. BEWICK, J. SHANG, AND W. F. FAGAN
[23] J. S. McLachlan, J. J. Hellmann, and M. W. Schwartz, A framework for debate of assisted migration in an era of climate change, Conservation Biology, 21 (2007), pp. 297–302. [24] L. W. Morrison, M. D. Korzukhin, and S. D. Porter, Predicted range expansion of the invasive fire ant, Solenopsis invicta, in the eastern United States based on the VEMAP global warming scenario, Diversity Distributions, 11 (2005), pp. 199–204. [25] J. Ni, A simulation of biomes on the Tibetan Plateau and their responses to global climate change, Mt. Res. Dev., 20 (2000), pp. 80–89. [26] S. Normand, R. E. Ricklefs, F. Skov, J. Bladt, O. Tackenberg, and J.-C. Svenning, Postglacial migration supplements climate in determining plant species ranges in Europe, Proc. Roy. Soc. B Biological Sci., 278 (2011), pp. 3644–3653. [27] C. V. Pao, Nonlinear Parabolic and Elliptic Equations, Plenum Press, New York, 1992. [28] A. Paradis, J. Elkinton, K. Hayhoe, and J. Buonaccorsi, Role of winter temperatures and climate change on the survival and future range expansion of the hemlock wooly adelgid (Adelges tsugae) in eastern North America, Mitigation Adaptation Strategies Global Change, 13 (2008), pp. 541–554. [29] F. Parker-Allie, C. F. Musil, and W. Thuiller, Effects of climate warming on the distributions of invasive Eurasion annual grasses: A South African perspective, Climate Change, 94 (2009), pp. 87–103. [30] C. Parmesan, Ecological and evolutionary responses to recent climate change, Ann. Rev. Ecology Evolution Systematics, 37 (2006), pp. 637–669. [31] C. L. Parr, E. F. Gray, and W. J. Bond, Cascading biodiversity and functional consequences of a global change-induced biome switch, Diversity Distributions, 18 (2012), pp. 493–503. [32] J. J. Polovina, J. P. Dunne, P. A. Woodworth, and E. A. Howell, Projected expansion of the subtropical biome and contraction of the temperate and equatorial upwelling biomes in the North Pacific under global warming, ICES J. Marine Sci., 68 (2011), pp. 986–995. [33] A. B. Potapov and M. A. Lewis, Climate and competition: The effect of moving range boundaries on habitat invisibility, Bull. Math. Biol., 66 (2004), pp. 975–1008. [34] F. J. Rahel and J. D. Olden, Asessing the effects of climate change on aquatic invasive species, Conservation Biology, 22 (2008), pp. 521–533. [35] M. Rejmanek, and D. M. Richardson, What Attributes Make Some Plant Species More Invasive? Ecology, 77 (1996), pp. 1655–1661. [36] A. Ricciardi and D. Simberloff, Assisted colonization is not a viable conservation strategy, Trends Ecology Evolution, 24 (2009), pp. 248–253. [37] M. Scheffer, M. Hirota, M. Holmgren, E. H. Van Nes, and F. S. Chapin, Thresholds for boreal biome transitions, Proc. Natl. Acad. Sci. USA, 109 (2012), pp. 21384–21389. [38] S. Scheiter and S. I. Higgins, Impacts of climate change on the vegetation of Africa: An adaptive dynamic vegetation modelling approach, Global Change Biology, 15 (2009), pp. 2224– 2246. [39] N. Shigesada, K. Kawasaki, and E. Teramoto, Traveling periodic waves in heterogeneous environments, Theoretical Population Biology, 30 (1986), pp. 143–160. [40] D. A. Smale and T. Wernberg, Extreme climatic event drives range contraction of a habitatforming species, Proc. Roy. Soc. B, 280 (2013). [41] R. W. Sutherst, Climate change and invasive species: A conceptual framework, in Invasive Species in a Changing World, H. A. Mooney, and R. J. Hobbs, eds., Island Press, Washington, DC, 2000, pp. 211–240. [42] J. C. Svenning and F. Skov, Ice age legacies in the geographical distribution of tree species richness in Europe, Global Ecology Biogeography, 16 (2007), pp. 234–245. [43] J. C. Svenning and F. Skov, Could the tree diversity pattern in Europe be generated by postglacial dispersal limitation? Ecology Letters, 10 (2007), pp. 453–460. [44] D. Y. P. Tng, B. P. Murphy, E. Weber, G. Sanders, G. J. Williamson, J. Kemp, and D. M. J. S. Bowman, Humid tropical rain forest has expanded into eucalypt forest and savanna over the last 50 years, Ecology Evolution, 2 (2012), pp. 34–45. [45] G. R. Walther, A. Roques, P. E. Hulme, M. T. Sykes, P. Pysek, I. Kuhn, M. Zobel, S. Bacher, Z. Botta-Dukat, H. Bugmann, B. Czucz, J. Dauber, T. Hickler, V. Jarosik, M. Kenis, S. Klotz, D. Minchin, M. Moora, W. Nentwig, J. Ott, V. Panov, B. Reineking, C. Robinet, V. Semenchenko, W. Solarz, W. Thuiller, M. Vila, K. Vohland, and J. Settele, Alien species in a warmer world: Risks and opportunities, Trends Ecology Evolution, 24 (2009), pp. 686–693. [46] X. Wang, On the Cauchy problem for reaction-diffusion equations, Trans. Amer. Math. Soc., 337 (1993), pp. 549–589.
SPREAD WITH A SHIFTING HABITAT EDGE
1417
[47] H. F. Weinberger, Long-time behavior of a class of biological models, SIAM J. Math. Anal., 13 (1982), pp. 353–396. [48] H. F. Weinberger, On spreading speeds and traveling waves for growth and migration models in a periodic habitat: Long-time behavior of a class of biological models, J. Math. Biol., 45 (2002), pp. 511–548. [49] H. F. Weinberger, M. A. Lewis, and B. Li, Analysis of linear determinacy for spread in cooperative models, J. Math. Biol., 45 (2002), pp. 183–218. [50] K. A. With, The landscape ecology of invasive spread, Conservation Biology, 16 (2002), pp. 1192–1203. [51] Y. Zhou and M. Kot, Discrete-time growth-dispersal models with shifting species ranges, Theoretical Ecology, 4 (2011), pp. 13–25.
