Spreading speeds of invasive species in a periodic patchy environment: effects of dispersal based on local information and gradient-based taxis

We consider a periodic environment with favorable and unfavorable patches alternately arranged in one dimension. In such heterogeneous environments, invasive animals may undergo not only random diffusion but also directed movement toward favorable patches. Here we propose reaction–diffusion–advection models in which the diffusion term is of either Fickian or Fokker–Planck type, the advection term is given by a gradient-based taxis, and the growth rate depends on both the environmental conditions and the population density. We first present hypotheses for the existence of a periodic traveling wave and a method to derive its spreading speed. Then this method is applied to a case in which the environment is homogeneous within each patch, so that taxis occurs at interfaces between different patches, causing accumulated distributions in favorable patches with density jumps at the interfaces. We show how the Fickian or Fokker–Planck diffusion and the gradient-based taxis interplay to determine the spreading speed.


Introduction
Environments in nature are generally heterogeneous as the consequence of natural processes or human activities causing habitat fragmentation. In such heterogeneous environments, the dispersal of organisms could be influenced by various factors such as spatially varying diffusion and directed movement toward more favorable habitat (taxis). At the same time, the rate of population growth would also change depending on the favorability of environments. In this article, we introduce reaction-diffusionadvection equations incorporating the above spatially dependent dispersal and growth processes, and investigate how these factors interplay in determining the spatiotemporal distribution of invasive species and its rate of spread. As for the spatial heterogeneity, we consider a simple periodically varying environment, in which a favorable patch and an unfavorable patch are alternately arranged in one dimension [36].
Many organisms from bacteria to mammals have abilities to migrate in response to stimuli or signals indicating foods, favorable or unfavorable habitats, prey or predators, etc., through various senses such as sight, hearing, smell, touch and so on. Since the pioneering work by Keller and Segel [17,18] on bacterial chemotaxis, various mathematical models in the framework of reaction-diffusion-advection equations have been proposed and used to investigate the taxis-induced pattern formations widely observed in not only microorganisms but also higher animals [1,2,6,8,9,13,14,23,28,29,31,34,38,41]. However, there have been relatively few works on the effect of taxis on the range expansion pattern of invading species and its spreading speed in heterogeneous environments, but see Kawasaki et al. [15], Maciel and Lutscher [24] and Berestycki et al. [3].
In the present work, we employ two types of reaction-diffusion-advection equations in one dimension in which diffusion terms are of the Fickian and Fokker-Planck types, and advection term is given by a gradient-based taxis: Fickian reaction-diffusion-advection equation (Fickian RDA) Fokker-Planck reaction-diffusion-advection equation (Fokker-Planck RDA) where n(x, t) is the population density, D(x) is the diffusion coefficient, u(x) is the advection velocity which is assumed to be proportional to the gradient of the favorability of environment f (x), i.e., u(x) ∝ d f (x)/dx, and R(x, n) is the per capita growth rate. We also assume that D(x), u(x) and R(x, n) are periodic with period L in x.

Fickian diffusion vs Fokker-Planck diffusion and gradient-based taxis
The difference in movement behaviors between the two types of diffusions was previously discussed in the framework of a one-dimensional random walk model by Skellam [40] (see also Okubo and Levin [28], [12,21,41]). If the transition probabilities P i,i+1 and P i+1,i (the probabilities of moving from point i to point i + 1 and from i + 1 to i, respectively) depend only on the environmental conditions at the point of departure (i.e., local information), the diffusion term is given by ∂ 2 D(x)n/∂ x 2 , which is designated here as Fokker-Planck diffusion after the terminology originally used by Turchin [41]. On the other hand, if both transition probabilities depend on the average condition between the two points (or the condition at the mid point) so that they are equal, the diffusion term is given by ∂{D(x)∂n/∂ x}/∂ x, which is designated as Fickian diffusion because this obeys the Fickian law. Since the Fokker-Planck diffusion term is rewritten as {D (x)n}, it can be regarded as the Fickian diffusion with an additional effective advection in the direction of decreasing D(x). Okubo and Levin [28] pointed out that the Fokker-Planck diffusion may be well suited to handle certain kinds of animal taxis (e.g., kinesis) and dispersal, while the Fickian type should be appropriate for the diffusion of physical properties and, in some cases, for biological diffusion in which the flux is always directed from high concentration to low concentration.
In the same framework of a one-dimensional random walk model as above, we derive the advection term induced by gradient-based taxis (Othmer and Steven [29]). Let us denote by f (i) the favorability at point i. The organism undergoes not only random diffusion (either Fickian or Fokker-Planck) but also directed movement at a rate proportional to the difference f (i +1)− f (i) in the direction where f (i +1)− f (i) is positive. Thus, the transition rates are assumed to be given, where k i,i+1 and k i+1,i are the transition rates due to random walk and the second terms represent the transition rates due to the directed movement. If k i,i+1 = k i+1,i = k, we have a diffusion advection equation where D = lim λ,τ →0 (λ 2 /τ )k and u(x) = lim λ,τ →0 (λ/τ )β{ f (i + 1) − f (i)}. λ is the distance between adjacent two points and τ is time interval of each step. When k i,i+1 and k i+1,i satisfy the condition that yields Fickian or Fokker-Planck diffusion, we have the Fickian or Fokker-Planck diffusion-advection equation, respectively.

Layout of the paper
In Sect. 2, it is shown that the Eqs. (1) and (2) can be transformed into an equation for which under some assumptions the existence of asymptotic spreading speeds, the existence of periodic traveling waves, and the existence of certain periodic traveling waves of the linearization of the differential equations and the relations between their speeds and the asymptotic spreading speeds are known. This means that accurate approximations to the spreading speeds can be found from tractable computations. The present work is devoted to carrying out such computations in order to obtain insight into the variation of the speed with the parameters of the problem. In Sect.

A model and its spreading speeds
We have presented two models (1) and (2) for the diffusion and advection of organisms. We shall make the following assumptions about the coefficients of these models. (1) and (2) have the following properties:

Hypotheses 2.1 The prescribed functions in
i. The function D(x) is positive, L-periodic, continuous and piecewise differentiable, and its derivative D (x) is uniformly bounded and piecewise continuous. ii. The function u(x) is L-periodic, integrable, uniformly bounded and piecewise continuous. iii. There is a positive constantn such that vi. There are an L-periodic uniformly bounded and uniformly positive function (x) and a positive constant γ such that (2).

Remark
In the numerical examples of this work, we shall set u(x) = αr (x), where r (x) := R(x, 0). In this case, the Hypotheses 2.1.ii becomes the requirement that r (x) should be L-periodic, continuous and piecewise differentiable, with its derivative r (x) uniformly bounded and piecewise continuous. In order to study (1) and (2) in a unified manner, we shall show how to reduce them to a simpler one by means of a trick.
is a uniformly positive uniformly bounded L-periodic solution of the equation where Because Eq. (2) can be written as which is of the form (1), a solution of (2) can also be put into the form (6) by replacing and making the change of variable (5). Of course, u(x) must also be replaced by u(x) − D (x) in (7).
This lemma will be proved in the Appendix 1.
In order to study the large-time behavior of a solution of the reaction-diffusionadvection Eq. (6) we shall apply some results of Weinberger [43]. We must, of course, make some assumptions on the prescribed functions.
We now introduce a discrete-time evolution model which is naturally related to the continuous-time model (6). For any function v 0 (x) let v(x, t) be the solution of (6) with the initial values v(x, 0) = v 0 (x), and define the time-1 operator Q by the condition that Because the Eq. (6) has no explicit t-dependence, it is also true that for any Thus we see that if we restrict the values of t to integers, any solution v(x, t) of (6) is mapped into the sequence of function v n (x) := v(x, n), and this sequence satisfies Conversely, if we are given a solution of this recursion, we can recover a continuoustime solution of (6) by letting v(x, n + t) be the solution of (6) with the initial value v n (x) and with 0 ≤ t ≤ 1. Because t is uniformly bounded in this process, it turns out that the large-time asymptotics of a solution of (6) are the same as the large-n asymptotics of a solution of the recursion (8).
Our results will be based on the following Lemma, which permits us transfer the results proved in Weinberger [43] to solutions of the Eq. (6).  [43]. This lemma will be proved in the Appendix 2. Theorem 2.6 of Weinberger [43] concerns the speeds of periodic traveling waves (PTW) of the recursion (8). (Note that the periodic traveling wave has also been termed as "traveling periodic wave" [39] and "pulsating traveling front" [5].) A rightward PTW of speed c of the continuous-time Eq. (6) was defined by Shigesada et al. [39] and Berestycki et al. [5] to be a solution of (6) in which the function W (y, x) has the properties 1. For each y, the function W (y, x) is L-periodic in x. 2. For each x, the function W (y, x) is nonincreasing in y, and a. lim y→+∞ W (y, x) = 0, and b. lim y→−∞ [π 1 (x) − W (y, x)] = 0, where π 1 (x) is a uniformly positive Lperiodic equilibrium solution of (6).
A leftward PTW of (6) of speed c is a solution of the form v(x, t) = W (−x − ct, x) where W satisfies the same conditions. Weinberger [43] gave a natural extension to define a rightward or leftward PTW of a recursion (8) to be a solution of (8) of the form v n (x) := W (±x − nc, x), where W (y, x) has the properties listed above. When Q is the time-one map of a continuous process, it is easily seen that a PTW of the recursion is simply the sequence obtained by setting t = n in a PTW of the continuous process. Once There exists a uniformly positive L-periodic equilibrium solution π 1 (x) of (6).
There also exist two numbers with −c * − < c * + such that every solution v(x, t) of That is, c * + is the rightward asymptotic spreading speed of any of a large family of initial value problems for (6) and c * − is the leftward asymptotic spreading speed. 2. The rightward spreading speed c * + can be characterized as the smallest c for which a rightward PTW of speed c exists. Also c * − is the smallest speed of a leftward PTW of (6). 3. c * + can also be characterized as the smallest speed c for which the linearization of (6) has a solution of the form with s > 0 and ψ(x) L-periodic. c * − is characterized as the smallest c for which (9) has a solution of the form w(x, t) = e −s[−x−ct] ψ(x) with s > 0 and ψ(x) L-periodic.
The statements of this Theorem can also be obtained by using the methods of Berestycki et al. [3][4][5].
As we shall see, Statement 3 of the Theorem 2.1 greatly simplifies the computation of the spreading speeds.

Fickian RDA equations for a periodic patchy environment
We consider a periodic patchy environment consisting of favorable and unfavorable patches with size l 1 and l 2 , respectively, and propose the Fickian RDA equation with a logistic growth function as below: x m indicate the left boundaries of the favorable patch located between m L and m L +l 1 . The diffusion coefficient D(x) and the intrinsic growth rate r (x) are d 1 (> 0) and r 1 (> 0) in the favorable patch and d 2 (> 0) and r 2 (< r 1 ) in the unfavorable patch, respectively. u(x) represents the advection velocity at position x, which means that organisms move right or left at speed |u(x)| when u(x) is positive or negative, respectively. In general, the functional form of u(x) may change with species depending on by what means and how far they sense the favorability of environments. As mentioned in the introduction, here we focus on the gradient-based taxis as given by where f (x) is the favorability of an environment at x and α the taxis sensitivity. As a simple and mathematically tractable candidate for f (x), we adopt the intrinsic growth rate r (x) := R(x, 0) [2,9,37]. Thus the advection velocity function is written as Because the functions D (x), u(x) and r (x) do not satisfy the assumptions of Hypotheses 2.1.i, ii and Remark below Hypotheses 2.1, maximum principle methods cannot be used to define the spreading speed and PTW. Instead, we define the spreading speed and PTW of the problem (10) as the limits as h → 0 of the the corresponding quantities for the one-parameter family In (12b) and (12c), we proposed a linear spline form forD(x) andr (x), which satisfy Hypotheses 2.1.i and Remark below Hypotheses 2.1. Since D(x), r (x) and u(x) are the limits as h → 0 ofD(x),r (x) andũ(x), respectively, we will obtain the solution of (10) as the limit of the solution of (12) when h → 0. As will be discussed in Sect. 5.1,D(x) andr (x) can be chosen from a more general class of functions.
In the examples of this and the next sections, the functions shows that the leftward spreading speed of such a problem is equal to the rightward spreading speed. Based on the above statement, we derive a rightward PTW solution and its spreading speed from the linearized equation (13) by using (14).
Let us introduce p(x) = e −sx g(x). Then (14) is rewritten as Substituting (15) into (13) yields By further introducing the variable q(x) =D(x) p −ũ(x) p , we have a system of first order differential equations, Because is L-periodic is rewritten in terms of p and q, Thus we only need to solve (16) with (17) for one spatial period. Here we focus on an interval, all parameters in (16) where 0 and l 1 are the left and right boundary points of the favorable patch, and 0 + and 0 − , or l + 1 and l − 1 , denote the limits as x approaches 0 or l 1 from right and left, respectively. Taking the consideration of The first equations in (19a) and (19b) indicate the ratios in population densities across the interfaces at x = 0 and x = l 1 , respectively. They are inversely related, and thus designated as κ F and 1/κ F , respectively. Since κ F > 1 always holds when α(r 1 − r 2 ) > 0, the gradient-based taxis causes a jump in density from unfavorable to favorable patches at the interface. The second equations in (19a) and (19b) represent the continuity in the flux at the interfaces. As special cases, (19a) is reduced to (20a) and Now we solve (16) on (0 − , L − ) by combining (17) for x = 0 − and the interface conditions (18a) and (18b). Then we have a dispersion relation between the two unknown variables, the speed of the PTW, c, at the limit of h → 0 and the damping coefficient, s, as below (see Appendix 3), where When there is a set of c > 0 and s > 0 that satisfy (21), c should be the speed of a PTW of (13) at the limit of h → 0. Theorem 2.1 states that the asymptotic speed c * is given by the formula c * = min s>0 c(s).
Then how can we derive the minimum speed? From (21) we have c = λ/s, which is rewritten as c = λ/Q(λ) by using (21a). Thus the speed of the PTW of (10), is given by We obtained the minimum value by using a numerical analysis software.

Numerical results
On the basis of the above theoretical study, we next examine how the spatio-temporal patterns of (10) and its spreading speed vary with the parameter values. We first carry out numerical simulations of equation (10) subject to the interface conditions, (19a) and (19b), with various initial distributions localized around the origin. In Fig. 1, the range expansions on the right half x axis after a sufficiently long time are shown, when α is 0, 0.5 and 2 while the other parameters are fixed as d 1 = 1, and dropping the primes for notational simplicity, the resultant dimensionless equation is given by (10) with d 1 = r 1 = μ = 1 [15]. Figure 1a shows the case without taxis (α = 0) in which the population density continuously changes at the interface (i.e., κ F = 1). In Fig. 1b, c where α = 0.5 and 2, respectively, the population density exhibits prominent jumps at interfaces (κ F = 2.5 for α = 0.5 and κ F = 38.3 for α = 2), while the population density within each favorable patch is gently convex or almost flat. Because a PTW of (12) is given by n : is L-preiodic (see (3)) and W (y, x) is nonincreasing in the first variable and L-periodic in the second one (see the statement below Lemma 2.2), a PTW n is nondecreasing in t for each x, and converges to b(x)π 1 (x) as t → ∞. Moreover, increasing t by t * := L/c * gives the graph with x increased by L. We see from Fig. 1 that the solution after a long time has the same properties as mentioned above. In particular, far behind the front, the graph is close to b(x)π 1 (x). Conversely, we can estimate the time t * when the graph is shifted by L, and then approximate c * from the formula c * = L/t * . From these and additional simulations, we also confirm that populations starting from various locally distributed propagules eventually either go to extinction or evolve to a unique PTW depending on parameter values. The condition for successful invasion (namely, n = 0 is unstable) is given by Q(0) > 0, where Q(λ) is defined by (21b, 21c). Figure 2 illustrates PTW speed c * as a function of the taxis sensitivity α for r 2 = 0.5, 0, −0.5, −1, −2, −3 while the other parameters are the same as in Fig. 1. The thick curves are the speeds of the Fickian RDA equation which are analytically derived from (22) with (21), and open circles are the speeds numerically calculated by L/t * . All curves are one-humped. Particularly, when r 2 is small, the speed c * sharply increases at first and then turns to decrease as α is increased. The large open circles indicated at α = 0, 0.5 and 2 for r 2 = −1 correspond to the range expansion patterns in Fig. 1a-c, respectively. In those three figures, the amplitude of variations in the population density between favorable and unfavorable patches increases as α increases. When α = 0.5, the population densities in the favorable patch is much higher than the corresponding density in the case of α = 0, and the population density in the unfavorable patch still remains at substantial levels though lower than that in the case of α = 0. These situations may contribute to the increased speed as seen in Fig. 2. However, when α = 2 in Fig. 1c, the population density in the favorable patch becomes close to the carrying capacity (r 1 /μ = 1), while that in the unfavorable patch approaches almost zero. In other words, organisms are trapped within the favorable patches so effectively that it would be harder for them to move to the next adjacent unfavorable patch, thereby leading to a decreased speed.
In order to see the effects of d 2 and l 2 in a broader range, we further obtain contour maps of c * on (d 2 , l 2 ) plane when α = 0, 0.5, and 2, while other parameters are kept the same as in Fig. 1. Figure 3a shows the results for the case in the absence of taxis (i.e., α = 0; see Shigesada et al. [39] for detail). The thick dotted curve indicates the boundary for population persistence, above which invasion fails. The boundary curve has both horizontal and vertical asymptotes which are given by lines l 2 = 1.09(≡ l c 2 ) and d 2 = 0.298, respectively (see dotted lines). When l 2 > l c 2 , c * shows a one-humped curve as d 2 increases and becomes zero when d 2 reaches a point on the boundary curve. As α increases from zero, the boundary curve shifts rightward, and thus the area for persistence is enlarged. For example, when α = 0.5, the vertical asymptote moves to d 2 = 1.51, while the horizontal asymptote is kept the same as l 2 = l c 2 = 1.09 (see

Fokker-Planck RDA equation 4.1 Formula for the PTW speed of the Fokker-Planck RDA equation
Here we consider the Fokker-Planck RDA equation with the logistic growth function as where D(x), r (x) are defined by (10b) and u(x) by (11). In a way similar to the case of the Fickian RDA equation described in Sect. 3, we consider a modification of (23) in which D(x), r (x) and u(x) are replaced byD(x),r (x) andũ(x) as defined in (12b), (12c) and (12d), respectively, As noted in Sect. 1.1 and Lemma 2.1, the above equation can be converted into a Fickian RDA equation, where the effective advection velocity due to the Fokker-Plank diffusion, −D (x), is newly added to the gradient-based taxis velocity,ũ(x). Thus we can apply the same analyses as done for the Fickian RDA equation (12) to (25) in which the advection velocityũ(x) originally present in (24) is replaced by −D (x) +ũ(x). As a result, we have the interface conditions corresponding to (19a) and (19b) in the Fickian RDA equation (see Appendix 4), and where κ r represents the density jump at the interface in the Fokker-Planck RDA model, which is equal to (d 2 /d 1 )κ F . (d 2 /d 1 ) comes from the effective advection caused by the Fokker-Planck diffusion based on local environmental information (see Appendix 4). Although κ F > 1 always holds, κ r could be either larger or smaller than one depending on the value of α(r 1 − r 2 )/(d 1 − d 2 ). Particularly, when α = 0 or r 1 = r 2 , we have which means that even in the absence of taxis, discontinuity in population density appears when d 1 is not equal to d 2 , in contrast to (20b) in the Fickian RDA model where density jump does not occur. The formula for the PTW speed of the Fokker-Planck RDA equation is given by (21) and (22) in which the equation for κ F in (21c) is replaced by κ r .

Numerical results
We first carry out numerical simulations of the Fokker-Planck RDA equation (23) with (26) and again confirm that initially localized populations eventually either go to extinction or evolve to a unique PTW depending on parameter values. The boundary curve for species persistence is given by Q(0) = 0, where Q(λ) is defined by (21b) and (21c) in which κ F is substituted by κ r . Figure 4a-c show the PTW solutions of the Fokker-Planck RDA equation with the same parameter sets as used in Fig. 1a-c, respectively. In Fig. 4a where α = 0, there appears density jump at the interface with ratio κ r = d 2 /d 1 = 1.2 so that the amplitude of variations in the population density between the favorable and unfavorable patches is larger than in Fig. 1a where the density changes continuously at the interface. When α = 0.5 and 2, the density jump is increased to κ r = 3.0 and 46.0, which are 1.2 times larger than the corresponding κ F s, though their spatio-temporal patterns are qualitatively similar to those in Fig. 1b, c, respectively. In Fig. 2, speed c * as a function of α for the Fokker-Planck RDA model is shown by the thin curves with closed circles. Comparing these thin curves with the thick ones for the Fickian RDA model, we can see that, for each r 2 , the thin curve is higher than the thick one at α = 0, whereas this relation is reversed when α exceeds around 0.5. This may be explained in such a way that since κ r is d 2 /d 1 (= 1.2) times higher than κ F , organisms with the Fokker-Planck diffusion are more strongly attracted to the favorable patches than those with the Fickian diffusion at sufficiently large α, resulting in decelerated speeds compared with the speeds of the Fickian RDA model. Figure 5a-c show contour maps of c * on (d 2 , l 2 ) plane, in which all parameters are the same as in the corresponding Fickian RDA model used in Fig. 3a-c. The closed circles in a-c correspond to the speeds of the PTWs shown in Fig. 4a-c, respectively. When α = 0 (Fig. 5a), the boundary curve for persistence appears in the upper left corner (thick dotted line). Note that, on the thin vertical lines where d 2 = 1 (i.e. d 1 = d 2 = 1), the speeds exactly coincide between the Fokker-Planck RDA and Fickian RDA models, because κ F = κ r = 1. In the region where d 2 < 1 and hence κ F > κ r , the speed c * of the Fokker-Planck RDA model (Fig. 5a) is lower than that of the Fickian RDA model (see Fig. 3a), and vice versa where d 2 > 1. Consequently, the area for population persistence in the region where d 2 < 1 is smaller in the Fokker-Planck RDA model than in the Fickian RDA model, and vice versa where d 2 > 1. As α increases from 0, the boundary curve for persistence (thick dotted line) approaches the y-axis, and vanishes when α ≤ 0.35 so that the species becomes able to survive for any positive set of (d 2 , l 2 ). On the other hand, the area for non-persistence in the Fickian RDA model moves to the right as α increases, as mentioned before (Fig. 3b,  c). Moreover, the above inequality relations in speed c * between the two models when α = 0 become not necessarily true. For example, in Fig. 5c (i.e., at α = 2), speed c * of the Fokker-Planck RDA model is higher than that of the Fickian RDA model (Fig. 3c) when d 2 < 1, and vice versa if otherwise. Notably, we found that, in the Fokker-Planck RDA model when α = 0 and d 1 < d 2 , there appears an unusual pattern such that the population density is higher in the unfavorable patch than in the favorable patch as shown in Fig. 6a, which corresponds to the open square in Fig. 5a. This is mathematically inevitable because the interface condition κ r = d 2 /d 1 = 0.8 is smaller than 1. However, when α = 0.5 and 2, κ F increases to 3 and 86.8 so that the density jump κ r = (d 2 /d 1 )κ F exceeds one to reach κ r = 2.4 and 69.4, respectively. As the result the PTW patterns become qualitatively similar to those in Fig. 4b, c.
Recently, Maciel and Lutscher [24] have investigated the effect of individual movement response to habitat edges (i.e., the interface between favorable and unfavorable patches in the present model) on the population persistence and spatial spread by using a reaction-diffusion equation with D(x) and r (x) as defined in (10b). They adopted two types of boundary conditions entailing discontinuity in density at the interface, which were derived by taking a limit of a random walk with habitat preference in movement direction as below [30] where a and 1 − a represent the probabilities that an individual at the interface moves to patch 1 and patch 2, respectively. Note that both (28) and (29) apparently differ from the discontinuity conditions in either (19) or (26) Fig. 4. In a, the population density is higher in the unfavorable patch than in the unfavorable patch. However, as the taxis sensitivity α is increased to 0.5 and 2, PTWs exhibit patterns qualitatively similar to those in Fig. 4b, c. The values of t, t * and c * are a 175.6 1.49 and 1.01, b 124.5, 1.07 and 1.40, and c 257.6, 2.28 and 0.66, respectively derived on the basis of mechanistic rule about individual movement behavior at x = 0, whereas the latter are deduced from the diffusion-reaction model with taxis which is induced by the spatial gradient in the growth rate. However, when there is neither gradient-based taxis (i.e., α = 0) nor biased movement (i.e., a = 0.5) at the interface, the interface condition in the Fokker-Planck RDA equation as given by (27) exactly coincides with (29) for a = 0.5, and hence speeds c * of the two models should be equal to each other when all other parameters are the same. They further examined the effect of a on the spreading speed c * , when patch preference a is a function of l 2 . Under these conditions, we cannot compare their models with ours as such. Nevertheless, their results seem qualitatively similar to those obtained from the Fokker-Planck RDA equation, but not the Fickian RDA equation, in the present study.

Discussion
We have previously studied a Fickian reaction-diffusion equation (i.e., (10) with u(x) = 0) to investigate the range expansion of invasive species in a periodic patchy environment [19,39]. There we numerically found a PTW solution and derived a heuristic formula for the speed. Subsequently Weinberger [43] rigorously proved the existence of the PTW for more general class of reaction-diffusion equations (see also, Berestycki et al. [3][4][5], Liang and Zhao [22]). In this article we extended our previous model to include an advection term while the diffusion term is chosen to be either the Fokker-Planck type or the Fickian type. In Sect. 2, Theorems for the existence of the PTW for general class of reaction-diffusion-advection equations are presented. It is also shown that when the growth function is of logistic type, the method to derive the PTW speed of the RD equation is applicable as such to the extended RDA equations. When the environment is homogeneous in each patch, the taxis occurs only at the interface of favorable and unfavorable patches, so that discontinuity in population density at the interface crucially influences the PTW speed. More specifically, the taxis acts to increase the PTW speed when the taxis sensitivity is moderate, but it turns to decelerate the speed when the taxis sensitivity becomes too large.

Interface condition
We assumed in Sects. 3 and 4 that D(x), r (x) and u(x) are the limits as h → 0 ofD(x), r (x) andũ(x)(= −dr (x)/dx) as given by the linear spline functions, (12b), (12c) and (12d), and obtained the interface conditions (19) and (26) for the Fickian RDA and the Fokker-Planck RDA models, respectively. However, it should be noted thatD(x) and r (x) could take various other functional forms within the region −h < x < h under the conditions that satisfy the Hypotheses 2.1.i, ii and Remark below Hypotheses 2.1. In such cases, the interface condition can generally be obtained by using the following equations (see Appendixes 3 and 4, and also Ovaskainen and Cornell [30]) The transitions (12b) and (12c) from D(x) toD(x) and r (x) tor (x) for −h ≤ x ≤ h can be written as the mollifications where ψ(x) = 1/2 for |x| ≤ 1, 0 for |x| > 1.
It is natural to look at the more general problem in which the mollifier uses a more general, and possibly smoother, function ψ(x) which satisfies the conditions ψ(x) ≥ 0, ψ(x) = 0 for |x| > 1, and D(x) andr (x) as defined by (31) with (32) are rewritten as: where Then the gradient-based taxis velocity is given bỹ Substituting (33a) and (34) into (30a), we obtain the interface condition at x = 0: Similarly, we use (30b) to have

Effective advection in the Fokker-Plank reaction-diffusion equation
We have investigated the effects of the gradient-based taxis, by which organisms are driven towards patch with larger r (x) at the interface. As noted in Sects. 1.1 and 4, however, even organisms without the gradient-based taxis can show an effective advection with velocity, −d D(x)/dx, if they undergo the Fokker-Planck diffusion, which depends only on local information. Thus let us consider the case that D(x) is a decreasing function of r (x), for instance, in analogy to the receptor-binding model used for bacterial chemotaxis [13], where γ is sufficiently small to make D(x) be uniformly positive.
Here we again assume that r (x) is given by (10b), which is the limit of h → 0 of r (x) defined in (12c). Then Fokker-Planck reaction-diffusion equation without taxis term is given by whereD(x) = D 0 /(1 + γr (x)) and the effective advection velocity is −dD(x)/dx = γ (1 + γr (x)) 2 dr (x)/dx. By substituting 0 toũ(x) in (30b), we have the interface condition in density as which becomes larger as either r 1 increases or r 2 decreases. Thus we could obtain accumulated distributions in the favorable patches with discontinuity jumps at the interfaces, despite the absence of an explicit advection term in Fokker-Planck RDA. However, the above model would not be suitable for macroscopic organisms which are capable of probing their surrounding environments and undergo directional movement.

Remaining problems for future studies
The gradient-taxis model presented in this article is concerned with the case that organisms have a relatively short sensing range at their current position. On the other hand, a higher organism with much longer sensing range may be able to integrate the favorability of environments within its sensing range. Thus we need to develop a model that addresses such long-range taxis.
Another limitation of the present models is that they do not take into account the population pressure which comes from density-dependent diffusion due to repulsion among individuals, which has been known to occur in various kinds of organisms [27,28,35]. Since the population density tends to be elevated by the presence of taxis as seen in Figs. 1b, c and 2b, c, the population pressure would have important effects. Thus we have been trying to extend our model to incorporate the population pressure.
Other remaining problems include extensions of the reaction-diffusion-advection model in heterogeneous environments to a two-dimensional space, incorporation of the Allee effect [10,16,20,32,42], multiple-species interactions [7,11,25,26] and so on, toward further understanding of how taxis influences the spatio-temporal pattern and the spreading speed of the periodic traveling wave under more realistic ecological contexts.

The principal tool in verifying the other hypotheses is the well-known Comparison Principle
Comparison Principle If two functions v (1) (x, t) and v (2) (x, t) have the properties ∂v (1) ∂t and if there is a nonnegative constant M such that (1) − v (2) ] when v (1) ≥ v (2) , for all t ≥ 0.
This result follows immediately from applying the maximum principle for parabolic equations to the function e Mt−x 2 (2) (x, t)] (see, e.g., Theorems 7 and 10 of Chapter 3 of Protter and Weinberger [33]).
By definition, this states that if 0 ≤ u (1) That is, Q is order-preserving, which is the second hypothesis in Weinberger [43].
The lattice L of Hypothesis 2.1.iii of Weinberger [43] consists of the set of translations by integral multiples of L, and P is the interval [0, L).
To verify Hypothesis 2.1.iv of Weinberger [43], we first define the function π 0 (x) ≡ 0, and observe that it is an equilibrium of the Eq. (6).
We define the function z(x) where (x) is the function in Hypothesis 2.1.vi and b(x) is the function (3). This change of variable takes the hypothesis into the condition that z(x) is uniformly bounded, uniformly positive, and L-periodic, and that We also observe that by Hypothesis 2.1.v wherem where is positive and small, satisfies the inequality We choose so small that the last factor is negative for 0 ≤ t ≤ 1. The Comparison Principle then shows that the solution v(x, t) of (6) with v(x, 0) = z(x) is bounded below by e [γ /2]t z(x). In particular, The Comparison Principle shows that v(x, n) is nondecreasing in n. The Comparison Principle also shows that if is so small that z(x) ≤n/ min x {b(x)}, wheren is the constant in Hypothesis 2.1.iii, then v(x, n) ≤n/ min x {b(x)}, which is uniformly bounded. Standard results on solutions of parabolic equations show that solutions of (6) have smoothness properties including equicontinuity. Therefore, the bounded increasing sequence of L-periodic functions v(x, n) converges to a positive continous L-periodic function, which we denote by π 1 (x). It is an equilibrium of the recursion (8), and also of the Eq. (6).
The strong maximum principle shows that if the initial values of a solution v of (6) are nonnegative, then v(x, t) > 0 for t > 0. If the initial function is also L-periodic, then v(x, t) is L-periodic in x. Therefore, v(x, 1) is bounded below by a multiple of z(x) with > 0, so that if π 0 = 0 ≤ v(x, 0) ≤ π 1 (x) then v(x, t) must converge to π 1 (x). Thus Hypothesis iv of Weinberger [43] is valid.
The standard results we have mentioned above also give the remaining two Hypotheses 2.1, so these Hypotheses follow from our Hypotheses 2.1.
Corollary 2.1 of Weinberger [43] requires additional conditions on the time-1 map M of the linearization of (6). Two conditions are required. Because Hypothesis 2.1.v makes the Comparison Principle states that if the initial values of a solution v of (6) and w of (38) are both v 0 (x), then v(x, t) ≤ w(x, t).
Setting t = 1 gives the condition This is one of the conditions of Corollary 2. we see that q(x, t) ≤ v(x, t) ≤ w(x, t) ≤ β for 0 ≤ t ≤ 1.
We use this inequality in the termmq in the equation for q to see that ∂q ∂ x − [S(x, 0) −mβ ]q ≥ 0 for 0 ≤ t ≤ 1.
We set t = 1 to see that Since e −mβ ≥ 1 − δ when is sufficiently small, this completes the proof of Lemma 2.2.
Thus, the following dispersion relation between speed c of a PTW and the damping coefficient s should hold where Multiplying the both sides of (39a) by exp − Then integrating the both sides of (44) over (−h, h) leads to As for the interval, (l 1 − h, l 1 + h), the same procedure as for the interval (−h, h) is applied to give lim h→0 A (l 1 −h,l 1 +h) = 1/κ F 0 0 1 .

Appendix 4: Derivation of (26)
Consider the linearized equation of (25), whereD(x),r (x) andũ(x) are defined by (12b), (12c) and (12d), respectively. Since the advection velocity, −D (x) +ũ(x), is uniformly bounded, all properties in Hypotheses 2.1 are satisfied if the equilibrium stateñ = 0 is unstable. Thus we can apply the same method as used in Sect. 3.1 to (50) to have the interface conditions and the speed of the PTW. Following (15), we putñ(x, t) = e −s(x−ct) g(x) = e sct p(x) as a PTW solution of (50). Then we have the same equations as (16) and (17), in whichũ(x) is replaced by −D (x) +ũ(x). Thus, substituting −D (x) +ũ(x) intoũ(x) in (46) of Appendix 3, we obtain By using n(h) = e sct p(h) in (15), we have the interface condition at x = 0 as the limit when h → 0, In a similar manner as in Appendix 3, the dispersion relation is given by (42) in which