EVOLUTION OF MIXED DISPERSAL IN PERIODIC ENVIRONMENTS

. Random dispersal describes the movement of organisms between adjacent spatial locations. However, the movement of some organisms such as seeds of plants can occur between non-adjacent spatial locations and is thus non-local. We propose to study a mixed dispersal strategy, which is a combina- tion of random dispersal and non-local dispersal. More speciﬁcally, we assume that a fraction of individuals in the population adopt random dispersal, while the remaining fraction assumes non-local dispersal. We investigate how such mixed dispersal aﬀects the invasion of a single species and also how mixed dis- persal strategy will evolve in spatially heterogeneous but temporally constant environment.


1.
Introduction. Dispersal of organisms has important consequence for population dynamics [3,14]. Individual organisms disperse to search for food, to breed, to avoid predators, and etc. Dispersal usually means the movement of organisms from one location to another, and there are many forms of movement, e.g., random dispersal and non-local dispersal. The underlying mathematical assumption for random dispersal is that organisms can only move to its immediate surrounding neighborhood and the transition probabilities in all directions are the same. One of the pioneering works in studying random dispersal is [38]. See [1,37] for more recent developments on the role of dispersal in biological invasions. Non-local dispersal, which is typical for the spatial spread of seeds of plants, assumes that organisms can travel for some distance and the transition probability from one location to another usually depends upon the distance the organisms traveled. For more background on non-local dispersal and related topics, see, e.g., [7,26,28]. There have been extensive studies on nonlocal diffusion models in recent years, almost exclusively for single species [2,4,5,6,8,9,10,12,22,24,25,29,31,32], but see [18,23] for two species. In this paper we propose to study a mixed dispersal strategy, which is a combination of both random dispersal and non-local dispersal. More specifically, we assume that a fraction of individuals in the population adopt random dispersal, while the remaining fraction assumes non-local dispersal. Our main goal is to investigate how the mixed dispersal affects the invasion of a single species and how the mixed dispersal strategies will evolve in spatially heterogeneous but temporally constant environment.
We assume that the environment will be spatially inhomogeneous and periodically varying in space. To facilitate our discussions, let p i > 0 for 1 ≤ i ≤ N . Set p = (p 1 , ..., p N ). We say that a function θ in R N is p-periodic if θ(x 1 , ..., x i + p i , ..., x N ) = θ(x) for every x ∈ R N and 1 ≤ i ≤ N . For the rest of the paper we will simply write θ(x + p) = θ(x) to indicate that θ is p-periodic in R N . Set D := (0, p 1 ) × (0, p 2 ) × · · · × (0, p N ). Define a nonlocal operator K : C(D) → C(D) by where w is p-periodic in R N , and the scalar function k(r) : [0, ∞) → [0, ∞) is assumed to be smooth, monotone decreasing and has compact support. Moreover, k(r) satisfies 1 ωN ∞ 0 k(r)r N −1 dr = 1, where ω N denotes the area of the unit sphere in R N . Note that k(|x − y|) is the probability of a non-local dispersing individual moving from location x to location y.
We first consider the invasion of a single species which adopts a combination of random dispersal and nonlocal dispersal, in an environment which is temporally constant and periodically varying in space. The formal derivations from Section 2 yield the following integro-partial differential equation model for a single species: where w(x, t) denotes the density of species at location x and time t, is the Laplace operator in R N which accounts for random dispersal of species, d is a positive constant which measures the total number of dispersal individuals per unit time, and the constant τ measures the fraction of individuals adopting random dispersal: τ = 0 corresponds to the scenario when all individuals adopt non-local dispersal, while τ = 1 represents the case when all individuals adopt random dispersal. We assume that τ ∈ (0, 1] unless otherwise specified. The function f represents the growth rate of the species and is assumed to be continuous differentiable in all components. We assume that the environment is periodically varying in space by imposing that f is p-periodic in x component, i.e. f (x + p, w) = f (x, w).
Consider the invasion of a single species governed by (1) subject to the periodic boundary condition: It is determined by the sign of the principal eigenvalue, denoted by λ 1 (d, τ, q), of the linear eigenvalue problem −d τ ∆ϕ + (1 − τ )Kϕ + qϕ = λϕ in R N , where function q(x) := −f w (x, 0) and is thus also p-periodic (see Proposition 3 for the existence and characterization of λ 1 (d, τ, q)). If λ 1 is positive, then the equilibrium solution w = 0 of (1)- (2) is stable and the species can not invade when rare. If λ 1 is negative, the equilibrium solution w = 0 is unstable and the species can invade when rare. Hence, the smaller λ 1 , the easier for the species to invade. Therefore, it is of interest to study the monotonicity of λ 1 with respect to parameters d and τ . Given any τ ∈ (0, 1] and non-constant q, one can show that λ 1 (d, τ, q) is strictly monotone increasing in d by the variational characterization of the principal eigenvalue. However, it is unclear whether λ 1 is also monotone in τ . In this connection, we first establish some general comparison results for λ 1 in terms of both d and τ .
The following result is an immediate consequence of part (i) of Theorem 1.1 for the case d 1 = d 2 .

Corollary 1. For any given d and non
Corollary 1 implies that when the period is suitably small, the smaller τ is, the smaller λ 1 is, and thus the easier for the species to invade. Biologically, small period corresponds to the scenario when the habitat is fragmented as the environment (e.g, habitat quality) has larger spatial variation. As smaller τ corresponds to the situation when the species has more tendency to adopt nonlocal dispersal, a consequence of Corollary 1 is that if the habitat is fragmented, it is advantageous for the species to adopt non-local dispersal in order to invade successfully.
It is natural to inquire whether λ 1 is always a strictly increasing function of τ . Both analytical example in Subsection 3.2 and numerical examples in Section 5 suggest that if the period is suitably large, λ 1 is not necessarily a monotone increasing function of τ . This suggests that when the habitat is not fragmented and the environment has small spatial variation, adopting random dispersal could potentially be more advantageous for the invasion of a single species.
To further study the effect of mixed dispersal strategy on the invasion of species, we adopt an approach which is similar to that of Hastings [15]. We consider an integro-partial differential equation model for two competing species which are identical in their population dynamics except for their dispersal strategies: both species disperse by mixed dispersal which is a combination of random diffusion and nonlocal dispersal, but with different random dispersal and non-local dispersal rates. To be more specific, we consider the two species competition system complemented with periodic boundary conditions where functions u(t, x), v(t, x) are the densities of two species, d u , d v > 0 are their dispersal rates, respectively. Biologically, τ 1 and τ 2 account for the probability of individuals of species u and v to adopt random diffusion, respectively; 1 − τ 1 and 1 − τ 2 account for the probability of individuals of species u and v to adopt nonlocal dispersal, respectively. The function a(x) represents the habitat quality and we assume that a(x) is non-constant, positive and p-periodic in R N to reflect the spatial inhomogeneity of the environment. We assume that 0 < τ 1 , τ 2 ≤ 1, unless otherwise specified. We first review some previous works related with system (6). Some of previous works focus on Neumann boundary conditions instead of the periodic boundary condition (7).
To describe our results, set It is easy to check that d * ≥ d * for any τ 1 , τ 2 ∈ (0, 1]; Moreover, d * = d * = 1 if τ 1 = τ 2 . Our next result can be stated as follows. Theorem 1.2. Suppose that a(x) is positive, non-constant, continuous and pperiodic, d u , d v are positive constants, and τ 1 , is globally asymptotically stable among all positive initial data; if d u /d v > d * , then (0, v * ) is globally asymptotically stable among all positive initial data.
When τ 1 = τ 2 , as d * = d * = 1, Theorem 1.2 implies that (0, v * ) is globally asymptotically stable if d v < d u , and (u * , 0) is globally asymptotically stable if d u < d v . This is consistent with results from previous work on random diffusion [11] and non-local dispersal [18,20]; i.e., the slower diffuser drives the faster diffuser to extinction. However, Theorem 1.2 contains more information: for example, if we regard d u /d v as a bifurcation parameter and let it vary from zero to infinity, two semi-trivial steady states will always exchange their stability, which suggests that the system has a branch of coexistence steady states for an interval of values of d u /d v . Note that such interval, if exists, must be a subset of If τ 1 = τ 2 , such interval shrinks to a point, i.e., d u /d v = 1. It will be of interest to investigate how many times (u * , 0) and (0, v * ) will exchange their stability and to study the structure of positive steady states. One special but interesting case is particularly noteworthy: Suppose that a(x) is positive, non-constant, continuous and p-periodic.
Biologically Corollary 2 means that non-local dispersal is preferred over random dispersal in such scenario, provided that max 1≤i≤N p i ≤ √ 2π which means that spatial variation of the environment is suitably large. It will be of interest to understand the case max 1≤i≤N p i > √ 2π. Our numerical simulation results from Section 5 suggest that if the period is suitably large, d u = d v and 0 < τ 1 < τ 2 ≤ 1, then (0, v * ) can be globally asymptotically stable. This shows that the evolution of the mixed dispersal strategy could depend on whether the habitat is fragmented, which echoes the discussions following Corollary 1.
The rest of this paper is organized as follows. Section 2 is devoted to the derivation of the mixed dispersal model for single species. In Section 3 we present some principal eigenvalue theory for mixed dispersal operator with periodic boundary condition and establish Theorem 1.1. Theorem 1.2 is established in Section 4, where some basic properties of the solutions of (6)-(7) are also collected. Section 5 is devoted to discussions and numerical simulations of eigenvalue problem (3) and the dynamics of (6)-(7).

2.
Derivation of the continuous model for single species. In this section, we follow Hutson et al. [20] to derive a continuous model for a single species in R which assumes both local and non-local movement. Following Hutson et al., we start with a discrete time and discrete space model and then let the sizes of time and space tend to zero. Divide the real line into intervals each with length δ and discrete the time into steps of τ . Let u(i, t) denote the density of the species in the interval [iδ, (i + 1)δ] and time t. We assume that the individuals may adopt either local or nonlocal movement. More precisely, we assume that for any interval a fraction γ of individuals move to the left or right of the interval each with probability 1/2, and the other fraction (i.e., 1 − γ) of individuals can jump to any interval on the real line, where constant γ ∈ [0, 1]. The change of the total number of individuals in the interval [iδ, (i + 1)δ] between time t and t + τ , which is given by is determined by two components: (i) Local dispersal. Individuals in the intervals [(i − 1)δ, iδ] and [(i + 1)δ, (i + 2)δ] will move to the interval [iδ, (i + 1)δ] with probability γ/2, and the total number of these individuals are given by (γ/2)u(i−1, t)∆x and (γ/2)u(i+1, t)∆x, respectively; individuals in the interval [iδ, (i+1)δ] will move out to its neighboring intervals with probability γ, and the total number of these individuals is given by γu(i, t)∆x. The change of individuals in the interval [iδ, (i + 1)δ] between time t and t + τ due to local transport is given by (ii) Nonlocal dispersal. Following Hutson et al. (pp. 487-488, [20]), we assume that the total number of individuals departing [iδ, (i + 1)δ] and arriving at [jδ, (j + 1)δ] is proportional to the population size in the interval [iδ, (i+1)δ] which is u(i, t)δ, the size of [jδ, (j + 1)δ] which is δ and the amount of transit time which is τ . Let α(j, i) be the proportionality constant. Then, the number of individuals arriving at [iδ, (i + 1)δ] through nonlocal transport is given by and the number of individuals departing [iδ, (i + 1)δ] through nonlocal transport is given by Dividing the above equation by δτ , we have We shall assume throughout that the rate of transition between the various patches, α(z, y) only depends on the distance between patches, i.e., α(z, y) = α(|z − y|). Set k(x) = α(|x|)/ρ, where the dispersal rate ρ represents the total number of the dispersing organisms per unit time. Then Note that η and ρ have different units: η = Length 2 /time and ρ = 1/time. Following Hutson et al., we introduce another parameter and replace k(z) by where L is the spread which characterizes the non-local dispersal distance, then u will satisfy The unit of d is 1/time and τ is dimensionless. Then, If we add population dynamics to (9), we arrive at which yields the integro-partial differential equation model (1) in the case N = 1 for a single species.

CHIU-YEN KAO, YUAN LOU AND WENXIAN SHEN
3. Principal eigenvalue for mixed dispersal operators. In this section, we first present some principal eigenvalue theory for mixed dispersal operators in periodic environment. Subsection 3.1 is devoted the proof of Theorem 1.1. In Subsection 3.2 we present some analytical example which suggests the non-monotonicity of the principal eigenvalue λ 1 (d, τ, q) with respect to τ . Consider the eigenvalue problem in R N , subject to periodic boundary condition ϕ(x) = ϕ(x + p), where d > 0 and τ ∈ (0, 1] are constants, and function q ∈ C(R N ) and is p-periodic. Let X be as in (8) and The following proposition, which is not difficult to prove, is about the positivity and compactness of the semigroup generated by Proposition 1. For any given 0 < τ ≤ 1, d(τ ∆ + (1 − τ )K) − q(·)I generates a strongly positive compact analytic semigroup on X.
Proof. It suffices to prove that K a self-adjoint operator on L 2 per (D). It follows from a result in [23] (see the proof of [23, Theorem 2.6]). In the following, we provide a simpler proof. Set It then follows that Proof.
It then suffices to prove that Observe that there are q n ∈ X such that q n ∈ C N (R N , R), q n − q X → 0 as n → ∞, and the partial derivatives of q n (·) at x 0n up to order N −1 are zero, where x 0n ∈ R N is such that q n (x 0n ) = min x∈R N q n (x). By [35,Theorem B], λ 1 (d, 0, q n ) exists. By [36, Lemma 3.1], λ 1 (d, 0, q n ) → λ 0 (d, 0, q) as n → ∞.

EVOLUTION OF MIXED DISPERSAL 2057
3.1. Proof of Theorem 1.1. We first derive some inequalities for periodic functions.
Lemma 3.1. Suppose that θ is p-periodic and continuous. Then where the equality holds if and only if θ is a constant function.
Proof. For p-periodic function u, we have where the equality holds if and only if θ(· + z) = θ(·) for a.e. z ∈ R N . This together with the continuity of θ implies that θ ≡ const.
Lemma 3.2. For any continuous p-periodic function θ, and equality holds if and only if θ is constant.

CHIU-YEN KAO, YUAN LOU AND WENXIAN SHEN
Hence, with either equality holds if and only if a λ = b λ = 0 for every λ = (0, ..., 0), i.e. θ is constant. Since we see that (14) holds, and equality holds if and only if θ is constant.

3.2.
Non-monotonicity of λ 1 in τ . If the condition max 1≤i≤N p i < √ 2π is violated, the conclusions in Corollary 1 may not hold anymore. To this end, assume that N = 1. Set k δ (z) = k(|z|/δ)/δ. For p-periodic function u, we have which can be proved by expressing function u in terms of its Fourier series. Therefore, for large δ, the eigenvalue problem (11) with k being replaced by k δ can be approximated by the following eigenvalue problem: where ϕ is p-periodic. Let λ(d, τ, ) be the principal eigenvalue of (19) with q(x) = cos 2πx p . The following result shows that λ(1, τ, ) may not be monotonically increasing in τ when p > 2π. Lemma 3.4. Suppose that q(x) = cos 2πx p . If p > 2π, then for 0 < 1, λ(1, 1, ) < λ (1, 0, ).

Two species competition model with mixed dispersals in periodic environment.
This section is devoted to the study of the two species competition system (6)- (7). In Subsection 4.1 we collect some basic properties of competition models with mixed dispersals. Theorem 1.2 is established in Subsection 4.2.
Proof. It follows from the arguments in [23,Theorem 3.2].
(2) It can be proved by the similar arguments as in (1).
(2) It can be proved similarly.

4.2.
Proof of Theorem 1.2. We first establish the non-existence of positive steady state of (6)- (7) when In this connection, we argue by contradiction. Suppose that (u, v) is a positive steady state of (6)- (7), i.e., u, v > 0 satisfy with periodic boundary conditions u(x) = u(x + p) and v(x) = v(x + p). Hence, where q = −a + u + v. We claim that q is non-constant. If not, suppose that q = C for some constant C. Integrating the equation of u in D, we see that C D u = 0.
Since u > 0 in D, C = 0. Hence, u satisfies Multiplying the above equation by u and integrating in D, we have Note that D |∇u| 2 ≥ 0, with equality holds if and only if u is constant. By Lemma 3.1, with equality holds if and only of u is constant. Since τ 1 ∈ (0, 1], the only possibility is that u is a constant function. Similarly, we can show that v is a constant function. As −a+u+v = 0, it implies that a is also a constant function, which contradicts our assumption that a is non-constant. This contradiction shows that q is non-constant. For the case 0 < τ 1 ≤ τ 2 ≤ 1, d u /d v < d * is equivalent to (4) and d u /d v > d * is equivalent to (5), with d 1 , d 2 being replaced by d u , d v , respectively. By part (i) of Theorem 1.1, as q = −a + u + v is non-constant, we see that if (4) holds, then . However, these conclusions contradict (39). This shows that for the case 0 < τ 1 ≤ τ 2 ≤ 1, system (6)- (7) has no positive steady state if The case 0 < τ 2 ≤ τ 1 ≤ 1 is identical to the case 0 < τ 1 ≤ τ 2 ≤ 1, so we omit the proof.
Next we show that if d u /d v < d * , then (0, v * ) is linearly unstable. We argue by contradiction and suppose that (0, v * ) is not linearly unstable. Then λ 1 (d u , τ 1 , −a+ v * ) ≥ 0. By the equation of v * , we have λ 1 (d v , τ 2 , −a + v * ) = 0. Therefore, The global asymptotic stability of (u * , 0) (when d u /d v < d * ) follows from the nonexistence of positive steady state and linear instability of (0, v * ) (see Proposition 7). Similarly, the global asymptotic stability of (0, v * ) (when d u /d v > d * ) follows from the nonexistence of positive steady state and linear instability of (0, v * ). This completes the proof of Theorem 1.2.

5.
Numerical simulations and discussions. In the previous sections, we analyzed the property of the principle eigenvalue for mixed random and nonlocal dispersal operators and the local stability of the two semi-trivial steady states in different scenarios. In order to understand more about global dynamic behaviors of the solutions for general parameter settings, we perform a series of numerical simulations in one dimension.
We use simple finite difference method [13] to obtain the solution numerically. For simplicity, we choose D = (0, p) where p is the period and define an uniform grid of points x j = j · h where 0 ≤ j ≤ N and N = p h . We use N = 400 in our numerical simulations. The second-order accuracy approximation of the second derivative is

EVOLUTION OF MIXED DISPERSAL 2067
The kernel used in the numerical simulations is defined in the following way. Let k(·) ∈ C ∞ (R) be defined by where C > 0 is chosen such that R k(x)dx = 1. For given δ > 0, let The second order approximation of integration of kernel term where k ij , u j are k δ (y j − x i ) and u(y j ) and w j is the weight of trapezoidal rule. For the eigenvalue problem (3), the discretization leads to a discrete eigenvalue problem and the first eigenvalue can be easily computed via Arnoldi's method [27]. For the competition model (6), the discretization leads to a system of ordinary differential equation. We then integrate in time by using Matlab built-in function "solver" which was designed to solve system of ordinary differential equations. The equilibrium results shown in the following figures are obtained when the difference between the solutions of two successive iterations is less than = 1.e − 10 or the iteration count reaches 30000.
In the last subfigure of Figure 3, the maximum of u and v are plotted with respect to d u to show the transition of stable steady states. In this example, the transition happens once (around d u ≈ 0.367) in the region d * d v < d u < d * d v rapidly. Whether this is a general behavior requires further investigation.
In Figure 4, we show how the equilibrium state of u and v of the competition model (6) varies with respect to d u for a larger period p = 4π. The rest of parameters and initial conditions are the same as the previous example. We also observe that if d u /d v < d * = 17/73, then (u * , 0) is globally asymptotically stable; if d u /d v > d * = 9, then (0, v * ) is globally asymptotically stable. In the last subfigure of Figure 4, the maximum of u and v are plotted with respect to d u to show the transition of stable steady states. In this example, the transition also happens only once (around d u ≈ 0.0175) in the region d * d v < d u < d * d v rapidly and it is very close to d * d v .
Comparing the first subfigure in Figure 3 and the third subfigure In Figure 4, two simulations have the same parameters except p. We see that (u * , 0) is globally asymptotically stable with p = 1 while (0, v * ) is globally asymptotically stable with p = 4π. This indicates that Corollary 1.4 cannot apply to general p.