A Boundary Value Problem for Integrodifference Population Models with Cyclic Kernels

Dedicated to Chris Cosner in admiration for his steady and inspiring efforts to bring together theory and application in mathematical biology and ecology. Abstract. The population dynamics of species with separate growth and dispersal stages can be modeled by a discrete-time, continuous-space integrod-ifference equation. Many authors have considered the case where the model parameters remain fixed over time, however real environments are constantly in flux. We develop a framework for analyzing the population dynamics when the dispersal parameters change over time in a cyclic fashion. In particular, for the case of N cyclic dispersal kernels modeling movement in the presence of unidirectional flow, we derive a 2N th-order boundary value problem that can be used to study the linear stability of the associated integrodifference model. 1. Motivation and background. Humans depend on biodiversity for trade, sustenance , pharmaceuticals, and entertainment, not to mention the ethical imperative of wildlife protection [13]. However, Earth's biodiversity is under serious threat. The International Union for Conservation of Nature lists 21,285 species as threatened, which is likely an underestimate due to insufficient information [2]. In fact, extinction rates are projected to increase as the human population continues to grow and put increased demand on Earth's resources [2]. Many sources agree that the primary cause of threat is habitat loss and degradation, followed by direct exploitation and competition with invasive non-native species [2, 15]. How can we ensure that threatened populations are able to survive in a limited habitat, or that non-native species will not overwhelm native populations? Population models are invaluable to conservation ecology because they help us predict the likely impact of environmental change on a given population. In this paper, we consider a population model for species that have separate growth and dispersal stages. Many plants demonstrate this pattern of growth and dispersal, as well as a variety of arthropods (insects, arachnids, and crustaceans) [7, 5]. We restrict our attention to populations living in a bounded habitat surrounded by

1. Motivation and background.Humans depend on biodiversity for trade, sustenance, pharmaceuticals, and entertainment, not to mention the ethical imperative of wildlife protection [13].However, Earth's biodiversity is under serious threat.The International Union for Conservation of Nature lists 21,285 species as threatened, which is likely an underestimate due to insufficient information [2].In fact, extinction rates are projected to increase as the human population continues to grow and put increased demand on Earth's resources [2].Many sources agree that the primary cause of threat is habitat loss and degradation, followed by direct exploitation and competition with invasive non-native species [2,15].How can we ensure that threatened populations are able to survive in a limited habitat, or that non-native species will not overwhelm native populations?
Population models are invaluable to conservation ecology because they help us predict the likely impact of environmental change on a given population.In this paper, we consider a population model for species that have separate growth and dispersal stages.Many plants demonstrate this pattern of growth and dispersal, as well as a variety of arthropods (insects, arachnids, and crustaceans) [7,5].We restrict our attention to populations living in a bounded habitat surrounded by uninhabitable regions.In particular, we are interested in populations subjected to a unidirectional flow, such as species like the southern bull-kelp in New Zealand that live in a current or river environment [1].This unidirectional flow will tend to bias the dispersal of the population in one direction, which, along with the bounded domain, gives rise to the drift paradox -the phenomenon of populations surviving in a limited habitat when they are being continually washed out of it [11] and the associated critical domain size-the minimal habitat size necessary for the population to persist [9,11].We are particularly interested in whether a population with known (experimentally determined) growth and dispersal behavior will survive or go extinct if introduced to a certain habitat.
Previous research in this area has largely focused on time-independent models, which assume the growth and dispersal behavior remain constant over time [9,16,11].However, in reality such features are unlikely to be constant, but may have seasonal change, random variation, or other time-dependence.In this note we consider the situation where the kernels change over time to reflect seasonal or annual variation.In particular, we consider the persistence or extinction of a population under cyclic kernel models.We also briefly relate this to the case of random kernels as considered in [4,8].
We consider organisms with separate growth and dispersal stages.Dispersal is assumed to be continuous in space and occurring over a fixed time interval, while growth is independent of space but depends on the local population density.Denoting n t (x) as the population density at time step t, the general time-independent model for the population growth is where Ω is the habitable domain, f is the density-dependent growth function, and K(x, y) is the dispersal kernel, which represents the probability that an individual starting at y will settle at x during the dispersal period.We assume f ≥ 0 is a smooth, monotonically increasing function with f (0) = 0. We consider continuous populations and assume that the properties of K are such that the operator defined by ( 1) is a positive operator on C(Ω), the space of continuous functions (which will be true for the specific advective kernels we consider).This integrodifference model was first studied in an ecological setting by [9].However, there have since been many extensions, including applications to predator-prey models [12], fragmented habitats [16], and stream environments [11].
When we consider temporal variation, both the growth and dispersal could have time-dependence, which leads to the general model where K t (x, y) is the dispersal kernel and f t (x) is the growth function at the t th time step.We are particularly interested in the case where the kernels K t have predictable time-dependence, such as seasonal variation, but the same framework has also been used to study models with random variation in parameters associated with growth and dispersal [4,8].
To determine a population's ability to survive in a new habitat we consider the local linear stability of the trivial equilibrium solution of (1).Namely, consider the linearization of (1) defined by where L is the Fréchet derivative at 0 of the right-hand side of (1).Assuming L is a strongly positive, bounded, compact operator, the Krein-Rutman theorem implies that L has a principal eigenvalue λ 1 > 0 [10].If λ 1 < 1, then the zero equilibrium is linearly stable and small perturbations will tend back to 0 (extinction).On the other hand, if λ 1 > 1 then the zero equilibrium is unstable, providing favorable conditions for the population to survive (persistence) [4].For λ 1 > 1, we do not consider more detailed behavior such as convergence to a nontrivial equilibrium state.
To motivate our analysis of the asymmetric kernel case for advective flow, we first review the situation for a symmetric kernel on domain Ω = (0, L).Following [12], if the population moves under Brownian motion during the dispersal stage with diffusion coefficient D and settles with rate α, then the associated dispersal kernel is a Laplace kernel of the form where a = α/D.To determine the associated principal eigenvalue of L we follow [16] and consider the eigenvalue equation where R = f (0).Differentiating twice with respect to x we have This implies an eigenfunction φ satisfies the equation The boundary conditions can be derived from (4) and (5) which yields the following Sturm-Liouville equation for the eigendata: Analysis of this system will produce the principal eigenvalue λ 1 to determine extinction or survival.This framework can also be used to study the critical domain size by setting λ = 1 and determining the associated value of L for which (7) admits a nontrivial positive solution [16].
If the population is subjected to a unidirectional flow of velocity v, such as in a river or wind environment, then it was shown in [11] that the dispersal kernel satisfies the differential equation which, assuming the probability an individual settles at x approaches zero as |x| → ∞ and that K approaches the same value coming from the left or the right of the starting position y, leads to the asymmetric Laplace kernel given by where and Note that if v = 0, this kernel reduces to the symmetric Laplace kernel (4).
In [8], the authors consider a time-dependent version of (1) where the kernels alternate between two cases K 1 and K 2 , with linearized growth rates R 1 = f 1 (0) and R 2 = f 2 (0), which could model predictable seasonal or annual variation in flow and growth rates.In the linearization of this alternating model, the population density after a time step with kernel K 1 is given by and after the subsequent time step with kernel K 2 we have Combining these two steps, we have We can represent this by a single equation that models both steps as where and each time step in (12) represents two steps in the original system.Since this combined model gives the population density after every other time step, the longterm dynamics will be the same.Since the kernel in ( 12) is time-independent, we can determine its principal eigenvalue λ 1 .Furthermore, since (12) represents two time steps, we can think of λ 1 = √ λ 1 as an effective principal eigenvalue for one time step.In particular, if λ 1 > 1, then λ 1 > 1, indicating persistence, and if λ 1 < 1, then λ 1 < 1, indicating extinction.

2.
Principal eigenvalue for N -stage kernels.In [8], the authors analyze the principal eigenvalue of ( 12) by deriving a fourth-order boundary value problem similar to (7), only now in the context of asymmetric kernels (9).In particular, their boundary-value problem can be used to numerically approximate or analytically solve for the principal eigenvalue λ 1 of the associated alternating kernel model (12).
Here we generalize the work of [8] to the case of N alternating kernels.In particular, for arbitrary asymmetric Laplace kernels, we derive a 2N th -order, linear, constant-coefficient, ordinary differential equation and 2N boundary conditions that allow one to detemine the principal eigenvalue of the associated model.Thus we may determine the conditions for survival for a wide range of populations that disperse in multiple stages, such as a variety of woody plants that have been shown to disperse according to alternating supra-annual schedules [6].
Suppose we have a population that disperses in N alternating stages, with dispersal kernels K 1 , K 2 , . . ., K N and corresponding linearized growth rates R 1 , R 2 , . . ., R N .Then inductively we see that the population density (near zero) after the (t + N ) th stage is where we have let x = x N +1 for convenience.We can represent ( 14) as a single integrodifference equation of the form where is the combined dispersal kernel and each time step in (15) represents N steps in the original model.Since the principal eigenvalue λ 1 for (15) corresponds to N periods of growth, we can think of as an effective principal eigenvalue for one step (e.g., one season or year).
Consider the eigenvalue equation for the combined model (15), letting φ N denote the eigenfunction: If we define . . .
then, recalling the definition of R and K, we see that φ N (x) solves (17) if and only if Now assume that each dispersal kernel K i is an asymmetric Laplace kernel (9) with parameters α i , D i , and v i and corresponding a i and b i given by (10).Differentiating each of the equations in (18) twice and using (8) for each K i we obtain the following set of ordinary differential equations for φ i : . . .
We can obtain a scalar differential equation for the eigenfunction φ N using the following lemma.
we have where C(N, k) is the sum of all products of the coefficients c i,1 and c i,2 (1 ≤ i ≤ N ) such that the second indices sum to k and no two first indices are equal, and where terms with an odd number of coefficients being multiplied are negated.
Before proving the lemma, we give an example and a brief remark to help elucidate the meaning of C(N, k).
Example 1.When N = 3, we have the system of equations According to the lemma, the resulting equation we get from this system is 0 = φ For example, the coefficient consists of a sum of all products of coefficients c i,1 and c i,2 in (22) such that the second indices sum to 2 and no two first indices are the same.Notice also how all terms with an odd number of coefficients are negative.
be the set of N -tuples of elements from Z 3 = {0, 1, 2} and let s i be the i th entry of s ∈ S N .Furthermore, let Define the parity σ(s) of s ∈ S N to be the number of nonzero entries s i in s.Then, where c i,0 = 1 for all i.This provides a closed form for the coefficients C(N, k) as well as a potential way of computing them.Note that C(N, 0) = 1 (empty product) and C(N, k) = 0 for k < 0 or k > 2N (empty sum), which agrees with our previous definition since there are no ways to multiply the c i,1 and c i,2 so that the sum of the indices is less than zero or greater than 2N and no two first indices are the same.
Proof of Lemma 2.1.We will proceed by induction on N .For N = 1, our system is which is the desired form.Suppose the lemma holds for the case of N .In the N + 1 case, we have the system of differential equations: . . .
Since the first N equations satisfy the conditions of the lemma, by our induction hypothesis we have an equation ( 21) for φ N .Recalling C(N, 0) = 1 we can write this as Thus we can take 2N derivatives of (27) and substitute in for φ 27), we can replace all derivatives of φ N in the above sum with derivatives of φ N +1 , giving which, by distributing the sum and reindexing, can be written as Recalling that C(N, k) = 0 for k < 0 or k > 2N and C(N, 0) = 1, we may write this as 0 = C(N, 0)φ which, after folding the first three terms into the appropriate sums and recalling that C(N, k) = 0 for k < 0, may finally be written as To complete the proof, we show that To see this, first recall that no two first indices may be the same in any of the coefficients in a single term of C(N + 1, k), thus each of its terms must contain either c N +1,1 or c N +1,2 or no coefficient with first index N + 1.All possible terms that do not contain such a coefficient must be products of coefficients of the first N equations such that the second indices sum to k and terms with an even number of coefficients are negated, the sum of which is by definition C(N, k).Furthermore, any term in C(N + 1, k) with c N +1,1 in it will be c N +1,1 times a product of coefficients of the first N equations such that the second indices sum to k −1, guaranteeing that the total sum of the second indices is k.Then, since multiplying by c N +1,1 changes the number of coefficients in the product from odd to even or even to odd, the sum of all terms in . Similarly, any term in C(N +1, k) with c N +1,2 in it will be c N +1,2 times a product of coefficients of the first N equations such that the second indices sum to k − 2, and multiplying by c N +1,2 changes the parity, so the sum of all terms in is the sum of these three collections of terms, as stated in (29).Thus (28) becomes which satisfies (21) for the case N + 1, completing our proof by induction.
We can apply Lemma 2.1 to the system (20) for N alternating kernels, where and φ 0 = φ N (due to the structure of (20)), which makes the inhomogeneous equation ( 21) homogeneous.The resulting equation is a 2N th -order differential equation solely given in terms of the eigenfunction φ N .We remark that for N = 2, the resulting fourth-order equation agrees with the alternating kernel model in [8].
To complete the eigenvalue analysis we need the associated set of boundary conditions.If the domain is Ω = (0, L), then we can use the definition of K i (x, y) in ( 18) and their derivatives evaluated at 0 and L to derive the set of boundary conditions: for i = 1, 2, . . ., N .We can convert all of these boundary conditions into equations exclusively in terms of φ N using the following lemma.
Lemma 2.2.For the system of differential equations given in Lemma 2.1, if we have the boundary condition at x = x 0 : for i ∈ {1, 2, . . ., N }, where m i is a constant, then this leads to a boundary condition in terms of φ N given by 0 = where B mi (N, i, k) is the sum of all products of the coefficients c j,1 and c j,2 for j = i + 1, . . ., N and m i = m i,1 such that the second indices sum to k and no two first indices are equal, and where terms with an odd number of objects being multiplied are negated.
To clarify, we illustrate with an example: Example 2. For N = 3 there are three cases: i = 1, 2, 3.For i = 1, we get a fifth-order boundary condition (recall that we let m i = m i,1 for clarity in the statement of the lemma).For i = 2, we get a third-order boundary condition Effective principal eigenvalue λ 1 for the two-step (solid) and three-step (dashed) models as a function of domain length, where can be used to provide useful information about a given population (in this case the critical domain size).

2.1.
Combinatorics of the coefficients.The results of Lemmas 2.1 and 2.2 give rise to some interesting combinatorics.For instance, one may quickly notice from a few examples that the number of terms in C(N, k), the coefficient for φ , is the same as the number of terms in C(N, 2N − k), the coefficient of φ (k) N .The coefficients B mi (N, i, k) also have symmetry in their number of terms.In fact, the number of terms in C(N, k) is given by T (N, k), the trinomial coefficient that gives the coefficient of , 14].There are many interpretations of this sequence, including the number of lattice paths from (0, 0) to (N, k) using steps (1, 0), (1, 1), or (1, 2) and the number of ordered trees having N + 1 leaves at level three and N + k + 3 edges [14].An interpretation that allows us to connect T (N, k) to C(N, k) is the number of ways to put k indistinguishable balls in N distinguishable buckets, where no bucket is allowed to have more than two balls [3].This relates to our definition of C(N, k) given in (23), since here we are summing over N -tuples where each entry is in {0, 1, 2} and the entries sum to k.This corresponds to distributing k indistinguishable "1"s into each of the N entries, where no entry may have more than two "1"s.The connection to trinomial coefficients inspires another way to write C(N, k): where [x k ] indicates taking the coefficient of x k in the expanded polynomial.This is the case because when we take cross-terms, no two first indices will be the same, and since the second index matches the power of x, the sum of the second indices in a cross-term will always equal the power of x it multiplies.Since the c i,j are negated, all cross-terms with an even number of c i,j s will be positive and all cross-terms with an odd number will be negative.From this standpoint, it makes sense that the number of terms in C(N, k) is T (N, k), because if we replace every coefficient in the above product with 1, each cross-term contributing to the coefficient x k will contribute 1 to the sum, thus counting the number of terms.We can also see the relationship between C(N, k) and T (N, k) in their recurrences.As we derived in the proof of Lemma 2.1, the recurrence for C(N, k) is and for N = 1, we have so we originally have one coefficient in front of each φ n term.The recurrence for the trinomial coefficients is with T (1, 0) = T (1, 1) = T (1, 2) = 1 [14].Since we do not combine any of the coefficients in (36), we see that if T (N, k) counts the number of terms in C(N, k), then T (N + 1, k) will count the number of terms in C(N + 1, k).Since their initial values match up, this is indeed the case.We can see both from the expression of C(N, k) given in ( 23) and from the definition of T (N, k) that the total number of terms in all of the coefficients will equal Similarly, we see that the number of coefficients in B mi (N, i, k) is also a polynomial coefficient, T (N − i, k), which gives the coefficient of This inspires the representation of B mi (N, i, k) as From this and (33) we see that the total number of terms in all coefficients of either boundary condition of order 2(N The numbers T (N, k) have the same recurrence relation as T (N, k), but with initial values T (0, 2) = 0 and T (0, 1) = T (0, 0) = 1.These initial values correspond to the relationship that for N = i, we have where there is one term in the coefficients for φ i and φ i , but no coefficient for φ i .
The relationship to polynomial coefficients is an unexpected but delightful result of Lemmas 2.1 and 2.2.

2.2.
Connection to random kernel models.In [8] the authors consider the N = 2 case (alternating kernels) as well as the random case where the kernel parameters are chosen from a given distribution.In the random case there is no longer a method to obtain a time-independent operator and hence the eigenvalue analysis no longer applies, but building on the work of [4,8] there is an analogue of the principal eigenvalue in the random case, denoted by Λ which can be used to determine persistence in the same fashion (i.e., Λ = 1 is the threshold between extinction (Λ < 1) and persistence (Λ > 1)).We refer the reader to [4,8] for more details, but briefly consider the following situation which can be analyzed using our framework.
Suppose instead of repeating N kernels in a cyclic fashion, we choose one of the N kernels at random, with equal probability ("N -sided coin-flip").Since we expect to choose an equal number of each of the kernels in the limit as t goes to infinity, we might expect the long term behavior of the two models to be roughly the same.In fact, [8] show that for the N = 2 case, it appears (based on computer simulations) that Λ closely matches the principal eigenvalue λ 1 for the alternating kernel case.
Using our results for the N -step model and computer simulations to approximate Λ for the N -sided coin-flip model, it appears that this pattern holds for N = 3 as well (see Figure 4).In fact, for the N = 3 case the analytic data seems to align even better with the random data.We conjecture that Λ = λ 1 for the N -sided coin-flip model and the cyclic N -step model (assuming a fixed set of kernels K 1 , . . ., K n ).If true, this relationship would allow us to use λ 1 , which can be computed from the results of the previous section, to study persistence of populations that follow the random N -sided coin-flip model.3. Summary.In summary, we have developed a framework for studying integrodifference models with cyclically repeating kernels.In particular, we derived a 2N thorder boundary value problem with 2N boundary conditions that allow us to determine the principal eigenvalue for a model with N alternating asymmetric Laplace dispersal kernels.This generalizes the case of one kernel [16] and two kernels [8] to the N kernel case.We illustrated this with an example of computing the principal eigenvalue as a function of the domain length for the case of two and three alternating kernels (Figure 2).We also connected this to recent work concerning random kernels, exploring the relationship between the alternating N -step model and the N -sided "coin-flip" model, in which one of N Laplace kernels is chosen randomly at each time step.We conjecture that the effective eigenvalue Λ for the coin-flip model will equal the effective principal eigenvalue λ 1 for the alternating model.Although we chose to use asymmetric Laplace kernels inspired by population dynamics in streams, our approach can be applied to other settings where the kernel satisfies an associated differential equation (e.g., Gaussian kernels).
Integrodifference population models are of great importance in analyzing the dynamics of species that have separate growth and dispersal stages.Developing the theory of time-dependence is essential to create more realistic models, providing ecologists with greater flexibility to choose models that are more accurate for the specific populations they study.This research develops a new tool for studying timedependent models, with the ultimate goal of aiding the conservation of valuable species and ecosystems.

2N k=0 T
(N, k) = 3 N since each element in S N = N i=1 Z 3 will contribute one term to a coefficient somewhere in the resulting differential equation (note that the S N,k partition S N ).