Introduction

As spatial ecology has developed, a great variety of mathematical modeling approaches have been used to study questions at various levels of complexity. Integro-difference equations, which feature a continuous space but discrete time formulation of population dynamics, have proven especially useful for studying questions about population-level processes and species interactions. For example, integro-difference models have been used to predict changes in gene frequency [37,38,39,52,59], and characterize species' spatial dynamics [22,23,24,25,28,29,30,35,46]. Because integro-difference equations often admit traveling wave solutions of various kinds, a primary focus in many of these studies has been spatial spread (e.g., expansion of a population or a favorable allele). Examples include scenarios in which population fronts can expand spatially in an accelerating fashion [31] and cases where one or more species can (or cannot) outrun the pace of environmental change [26,34,61]. The reader is referred to the monograph by Lutscher [40] for a thorough review on integro-difference equations and applications.

Here, we adopt a very different perspective in that we use an integro-difference formulation to study nonspreading solutions. Roughly speaking, a nonspreading solution is a solution which persists with virtually bounded extent for all generations in an unbounded domain. Such a solution describes 'invasion pinning' that has been investigated for coupled ordinary differential systems in a discrete (patch) environment (see Keitt et al. [27] and references therein). Similar results can emerge for partial differential equations when the focus is on gap-crossing ability in heterogeneous landscapes, leading to 'geographic range margins' beyond which the species cannot spread [16]. Related results for gap crossing in integro-difference equations can be found in Musgrave et al. [43].

As we discuss below, the existence of nonspreading solutions in a homogeneous environment hinges on the presence of an Allee effect and overcompensation. An Allee effect arises when the per-capita birth rate increases as a function of population density when population density is small. Allee effects may occur via a great many biological mechanisms [1,2,7,8,9,12,13,14,15,21,36,42,49,53], and they have been studied in connection with integro-difference equations in the context of spatial spread [32,58]. A special kind of Allee effect, termed a strong Allee effect, occurs when there is a critical population density below which extinction occurs.

Mating failure, which can arise through mechanisms like pollen limitation and reproductive asynchrony [9,21], has been linked to strong Allee effects in diverse biological systems. For example, in evergreen bagworms (Thyridopteryx ephemeraeformis), the intensity of a strong Allee effect arising from mating failure is a function of climate, and this spatial variation leads to a hard geographic boundary for the species [41,50].

Overcompensation in population biology refers to phenomena in which density-dependent processes do not yield a smooth approach to carrying capacity, and overcrowding causes an overly dense population to decrease below carrying capacity, sometimes dramatically, rather than slowly declining to carrying capacity. Such imprecision in density dependence is often critical to the formation of cycles or chaotic dynamics in population models, and there is particular attention to the strength of overcompensation as a feature of the dynamics. One example of overcompensation in an ecological system is work by Symonides et al. [56] who demonstrated that overcompensation in germination success leads to cycles in the annual plant Erophila verna. Overcompensatory population crashes have been widely studied in small mammal species, where fast population growth rates, coupled with various combinations of parasitism, overexploitation of resources, and increased predation, are linked to the emergence of population cycles (e.g., [3,18]). Note that in studies of herbivory, overcompensation has a different definition that focuses on regrowth stimulated by herbivore feeding damage. That usage is not relevant here.

We consider the spatial-temporal dynamics of a population governed by the integro-difference equation

where u n (x) is the density of individuals at point x and time n, g(u) describes density dependent fecundity, and k(x-y) is the dispersal function, which depends upon the distance |y-x| between the location of birth y and the location of settlement x. Model (1) describes that individuals at location y generate g(u n (y)) offspring and then die and these offspring disperse to location x with the probability k(xy). We will assume that g(1) = 1, so that the population has an equilibrium at the carrying capacity u n (x) ≡ 1, and that k(x) decays at least exponentially fast near ±∞ so that the probability that an individual travels a very long distance is exponentially small.

In the case that small populations grow (i.e., g (0) > 1) and the reproduction function exhibits no Allee effect (i.e., g(u) ≤ g (0)u), with and without overcompensation, the population will spread at a constant asymptotic spreading speed that can be characterized as the slowest speed of a class of traveling waves (Weinberger [60], Li et al. [33]). In this case, the wave is pulled by the leading edge of invasion. When the reproduction function produces overcompensation, oscillations are generated in the population density behind the wave front (Bourgeois et al. [4,5],

Li et al. [33]). Constant spreading speed also occurs if g(u) exhibits a strong Allee effect (i.e., g (0) < 1) and g(u) is increasing (i.e., there is no overcompensation). In this case, the spreading speed is the unique speed of traveling waves connecting zero and the carrying capacity (Lui [39]), and the sign of the wave speed is the same as that of 1 0 [g(u)u]du (Wang et al [58]).

If 1 0 [g(u)u]du > 0, the traveling wave moves forward, if 1 0 [g(u)u]du < 0, the traveling wave moves backward, and if 1 0 [g(u)u]du = 0 the traveling wave is stationary. The wave speed depends on the forward pushing force developed by the high-density populations above the Allee threshold behind wave front as well as the backward pulling force generated by the lower-density populations below the Allee threshold along the leading edge of invasion.

Fluctuating invasion speeds can be generated by a strong Allee effect and strong overcompensation (Sullivan et al. [55]). Strong overcompensation in general produces large spatiotemporal variation in density behind the invasion front and thus, variation in the strength of the push, leading to oscillating spreading speeds. As pointed out in [55], where the population density is smaller than the Allee threshold along the leading edge of the invasion, the population declines before the next time step. Populations above the Allee threshold will grow until a maximum population is reached and overcompensation causes a reduction in growth. If overcompensation is strong enough, they will return from a high level to a low level resulting in cyclical variability in the pushing strength of the wave.

In this paper, we further study the effects of a combination of a strong Allee effect and strong overcompensation. We will show that such a combination can produce biologically meaningful nonspreading solutions that are robust in that they occur in solid regions of parameter space for (1). We will demonstrate the existence of nonspreading solutions with a variety of spatiotemporal patterns. One of our novel findings is the existence of nonspreading solutions that oscillate in both density and spatial range. Here in the long run, the oscillating forward pushing force developed by overcompensation is balanced by the backward pulling force from populations below the Allee threshold, leading to persisting nonspreading solutions. It should be pointed out that, as discussed above, for the case of no overcompensation, there exists a traveling wave with zero speed if 1 0 [g(u)u]du = 0; however this condition is not robust and is not satisfied with a slight change of model parameters. It is interesting to note that simple scalar integro-difference equations produce biologically meaningful nonspreading solutions. This is a phenomenon that scalar reaction-diffusion equations cannot produce. One consequence of the existence of nonspreading solutions concentrated on effectively bounded domains is the formation of multiple population 'patches' separated from each other in space. This will be also explored in the present paper.

Model Formulation

We will model growth using the two parameter function presented by Vortcamp et al. [57]. In this growth function, a represents the Allee threshold and r represents a parameter controlling the strength of overcompensation. By an appropriate scaling of u the carrying capacity can be assumed to be 1. The growth function used in model ( 1) is then

It can be shown that

The maximizer of g a,r (u) is given by

The resulting expression for g a,r (u max ) does not simplify into a compact form. By noting the signs of g a,r (a) and g a,r (1), we see a < u max < 1 if r > a 1-a .

Increasing r for fixed a increases the maximum value of g a,r (u) , increases g a,r (a) , and decreases g a,r (1). Conversely, increasing a for fixed r, decreases the maximum value of g a,r (u), decreases g a,r (a) , and increases g a,r (1).

In this parametrization, the shape of g a,r (u) is more sensitive to the parameter a than r, at least in range of parameters of interest to our study. As evidence for this, graphs of the partial derivatives of g a,r (u max ) with respect to a and r are shown in Fig. 1. g a,r (u max ) is the maximum value attained by g a,r (•). We see that ∂ ∂ a g a,r (u max ) is approximately two orders of magnitude greater than ∂ ∂ a g a,r (u max ). In Fig. 2 the effects on the graph of g a,r (u) of increasing a and r from a base case are shown. 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.5 1.0 1.5 2.0 2.5 a=0.4 , r=5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.5 1.0 1.5 2.0 2.5 a=0.5 , r=5 0.0 0.2 0.4 0.6 0.8 1.0 1.2 0.5 1.0 1.5 2.0 2.5 a=0.4 , r=6 Figure 2: The graph of g a,r (u) and y = u showing the effects of varying a and r.

Essential extinction occurs when severe overcompensation causes large populations to fall below the Allee threshold. Mathematically this is equivalent to image of the maximum value of g a,r (u)

being less than the Allee threshold [51]. The region in parameter space producing essential extinction is indicated in green in Fig. 13(a). As discussed in the introduction, for monotone growth functions with an Allee effect the sign of wave speed is equal to the sign of 1 0 [g a,r (u)u] du [58]. While this is not necessarily true for non-monotone functions, we can say the overall growth is weak if the integral is negative. The region in parameter space with weak growth is indicated in blue in Fig. 13(a).

It has long been known that the shape of the dispersal kernel, particularly its kurtosis, can have a profound influence on the spreading dynamics in an integro-difference equation [31]. To model dispersal with varying kurtosis, we will use the generalized Gaussian distribution [44] centered at the origin with standard deviation 1. Spatial coordinates can trivially be rescaled to satisfy that the standard deviation is 1 without altering the dynamics of the integro-difference equation.

The kurtosis of the distribution is controlled by the parameter η. While the distribution and all its moments are well defined if η > 0 , it is only exponentially bounded if η ≥ 1. The probability density function used in model (1) for dispersal is then

and Γ(•) refers to the gamma function.

Kurtosis, which is defined as the ratio of the fourth moment to the square of the second, gives a measure of the "fatness" of the tail of the distribution. Leptokurtic distributions (0 < η < 2) can be thought of as having most individuals disperse very small distances with a few individuals dispersing extreme distances in such a way the standard deviation remains fixed. Conversely, platykurtic distributions (η > 2) can be thought of as most individuals dispersing about the same distance. When η = 1, 2, and ∞ , the commonly used Laplace, Gaussian and Uniform distribution are recovered as is shown in Fig. 3. The spatial model is specified by model ( 1) with the definitions of g a,r (u) and k η (x) previously outlined. We model a unimodal symmetric initial conditions with a width and height parameter using

where p 0 is the maximum density and w 0 is the width of the support. The solution set {u n (x)} ∞ n=0 is thus fully specified by the five parameters a, r, η, p 0 , w 0 .

To numerically generate the solution set, we uniformly discretize space using a step size δ = 0.005 and use conv in Matlab to compute the accelerated convolution. We use the symmetry of u n (x) about x = 0 to further accelerate calculations. Both the vector representing k(x) and u n (x) are clipped where they fall below 10 -4 . The Matlab code can be viewed on https: //github.com/glotto01/theoretical-ecology.git. Numerical experimentation showed us that decreasing δ or the clipping threshold did not alter results. For example, Fig. 8 , was recreated using a clipping threshold of 10 -5 and δ = 0.0025 with identical results.

Nonspreading Solutions

In contrast to integro-difference equations with Allee or overcompensation effects considered separately [33,39,58], we are able to find solid regions of parameter space with solutions where the population persists but is effectively confined to a limited region of space. For example, in Fig. 4 we see a solution converging to a stable equilibrium where the population is effectively limited to -4 ≤ x ≤ 4. More complex behavior such as period-two and non-periodic nonspreading solutions can be observed as well (Fig. 5, 6). These unimodal nonspreading solutions can act as basis for solutions consisting of multiple patches, as is shown in Fig. 7. Throughout this paper, we define the spatial extent of generation t to be the distance from the left most point where u t (x) = a to the rightmost point where u t (x) = a.

In the following sections we will explore how nonspreading behavior depends on parameters and initial conditions, how two nonspreading solutions interact when superimposed, and finally how isolated patch like solutions can emerge from a noisy but nearly uniform initial density.

-4 -2 2 4 0.2 0.4 0.6 0.8 1.0

(a) density curves

5 10 15 20 25 30 2.0 2.2 2.4 2.6 2.8 generation spatial extent (b) spatial extent -4 -2 2 4 0.2 0.4 0.6 0.8 1.0 1.2

(a) density curves -4 -2 2 4 0.2 0.4 0.6 0.8 1.0 1.2

(a) density curves -10 -5 5 10 0.2 0.4 0.6 0.8 1.0 1.2 Figure 7: The population density curves for a solution forming "two" patches. Blue is the initial density, gray are the transients, and black is the apparent steady state. Parameters used are a = 0.61, r = 8, η = 5, p 0 = 0.8, w 0 = 25 3.1 Bifurcations for Nonspreading Solutions The Matlab functions used can be found in the one_parameter_orbit folder in https://github. com/glotto01/theoretical-ecology.git. In Fig. 8-12 we present orbit diagrams with respect to each of the parameters. With the exception of the parameter being varied, the other parameters are held at a = 0.61, r = 8, η = 5, p 0 = 1, w 0 = 6 . The bifurcations around this set of parameters are typical of those made for other choices based on our extensive simulations. The x-axis is the bifurcation parameter, and the y-axis is the spatial extent of the density curves u n (x) for 800 ≤ n ≤ 1000. 3000 uniformly spaced sample points in the bifurcation parameter are used to create the plot (Since η is half-log plot the actual sample points are geometrically spaced).

As can be seen in Fig. 456, typical time scales for transients are on the order of tens of generations but this can increase dramatically for parameters near bifurcation points (e.g. near a = 0.58 in Fig. 8). To insure the choice of 800-1000 was sufficient we computed the orbit diagram in Fig. 8 using 1600-1800 and found it to look identical to that presented here.

In Fig. 8 we see a period doubling bifurcation in the parameter a. For values of a between 0.603 and 0.64, we see a single period 1 solution emerge; for values between 0.58 and 0.603, a period-two solution emerges; and for values approximately between 0.575 and 0.58 higher order periodicities occur. Extinction occurs for small and large values of a, and it can be seen that regions of extinction are intermingled with nonspreading solutions for a between 0.58 and 0.6. While in the figure it appears that regions of extinction and survival overlap, that is an artifact of the point size used in plotting.

In Fig. 9 we show the orbit diagram for the parameter r. We see a period one nonspreading solution transition to a period-two, followed by extinction. In Fig. 10 we show the bifurcation behaviour for η, which is the parameter controlling the kurtosis of dispersal. We see that extinction occurs for leptokurtic dispersal (η < 2), period-two solutions occur for η slightly higher than 2 and less than 4, and a period 1 equilibrium for highly platykurtic dispersal when η > 4.

In the orbit diagrams for initial conditions, p 0 and w 0 (Fig. 11, 12), we see the solution is either attracted only to the period one orbit associated with that parameter set or to extinction. The basin of attraction for the period one orbit is fairly insensitive to w 0 , with widths ranging from 4 to 8 all producing the stable nonspreading solution. It should be noted that w 0 differs from spatial extent, as spatial extent measures the distance between where the population equals the Allee threshold whereas w 0 measures the support of the initial condition. 0.6 0.8 1. 2 4 6 8 10 20 40 60 80 100 0 1 2 3 4 η spatial extent Figure 10: Orbit diagram for parameter η with a = 0.61, r = 8, p 0 = 1, w 0 = 6. Note that the scale on the x-axis is logarithmic.

Two Parameter Bifurcations

The Matlab codes used in this section can be found in the folder 2parameter_bifurc on https: //github.com/glotto01/theoretical-ecology.git.

The nonspreading solutions exist along a fairly narrow band in parameter space centered around the curve where 1 0 [g a,r (u)u] du = 0 (see Fig. 13). As was shown in Fig. 1, the maximum value of g a,r (u) is more sensitive to the parameter a than r and the narrowness can be considered an artifact of this parametrization.

To identify regions in parameter space with different qualitative behavior, we divided the region of the ar plane depicted in Fig. 13(b) into a 100 × 100 grid. For the values of a, r on the grid, iterations are computed using η = 5, p 0 = 1, w 0 = 6. The following are used as criteria for classification:

• If for some n, the maximum value of u n (x) is less than a the solution is classified as extinction (gray).

• If periodicity is detected the solution is classified as nonspreading (blue). We will discuss how periodicity is detected below.

• If periodicity is not detected within 500 generations then it will be classified as spreading (yellow) if the ratio of the spatial extent for u 500 (x) to that of u 250 (x) is greater than 1.5, otherwise it will be classified as nonspreading. The justification for this threshold is discussed in Appendix A. To detect periodicity, the past values of the spatial extent are scanned for repeats. If a repeated value of spatial extent (to a tolerance of 10 -5 ) is detected, the density curve of the corresponding generation is compared to the present density curve. If the maximum absolute difference of the two density curves is less than 10 -5 it is classified as periodic.

Interaction of Two Patches

The Matlab codes used in this section can be found in the folder two_patch on https:// github.com/glotto01/theoretical-ecology.git.

In this section we will refer to a single unimodal equilibrium as a patch. We will use u p (x)

to refer to a single patch centered at the origin. We focus on the parameters a = 0.61, r = 8, η = 5 but the results for these parameters are qualitatively similar to that of other parameters possessing a single non-periodic patch solution based on our extensive simulations. We did not systematically investigate patch interaction for periodic solutions.

It's worth mentioning that for other parameters possessing a patch equilibrium, the shape, height, and width is similar to that in Fig. 4. Recall the dispersal kernel is scaled so the standard deviation is 1, this means it is impossible for the width of an equilibrium to be less than several units, due to taking a convolution with a probability function with a width of several units. The reason that equilibria with a width much larger than a few units do not occur, is because the growth parameters in the range that we are discussing exhibit essential extinction. Since the non-spatial model with essential extinction goes extinct almost surely [51], it is reasonable to assume that large spatially uniform densities would be unstable.

To investigate how two such patches interact with each other, we consider initial data in the

where d is the parameter controlling the separation of the patches. We find, to the limits we are able to effectively explore with a desktop simulation, that the population goes extinct for all values of d less than about 8. The time to extinction increases extremely rapidly with the parameter d as can be seen in Fig. 15. We terminated at d = 7.866 as the population did not go extinct in 500,000 generations and computational time became prohibitive. Here we are defining the time of extinction as the generation when the maximum value of the density falls below the Allee threshold.

We are not able to determine if mathematically stable equilibria exist for d > 8 or if they are just extraordinarily long lived transients. Since the dispersal kernel k η (x) falls off superexponentially (nonspreading solutions are not known to occur for η < 1), the overlap of the two patches presumably falls off super-exponentially in d. This would explain why the time to extinction increases so rapidly in d. Biologically the distinction between a true mathematical equilibrium, and an extremely long lived transient may not be as important of a distinction.

3 4 5 6 7 1 10 100 1000 10 4 10 5 d time to extinction

Patch Formation with a Stochastic Initial Condition

The Matlab codes used in this section can be found in the folder stochastic_initial on https: //github.com/glotto01/theoretical-ecology.git.

We wish to simulate a spatially stochastic initial condition with spatial correlation length of L scale . L scale can be considered as the length at which statistical correlations in density diminish to insignificant levels. The purpose of this is to examine the possibility of pattern emergence from a perturbed uniform initial density.

To accomplish this we generate N random numbers, y 1 , y 2 , • • • , y N , uniformly distributed on

N is chosen so that (N -1) L is the desired domain size. The initial density is then the linear interpolation of the points (i -1)L, y i for i = 1, 2, • • • , N.

We limit our study to the parameter values a = .61, r = 8, η = 5. The results are qualitatively similar to that of other parameters producing only a non-periodic nonspreading solution.

We did not systematically investigate the case for parameters producing periodic nonspreading solutions.

To study the effects of L scale on patch formation we:

• Used a domain of length 500 initiated as described above.

• Iterated until the maximum absolute difference of successive density curves is less 10 -5 .

• Counted the number of patches formed by integrating the total population and dividing by the population of a single patch. The population of a single unimodal patch is u p (x)dx = 3.695.

• For values of L scale = 0.01, 0.05, • • • , 50, 100 , twenty trials were completed for each value.

• 90% confidence intervals are computed using the assumption of normality (Student T distribution). It should be noted that the sample standard deviation for patch formation was about 3 patches independent of L scale .

The length scale of a single patch is ∼ 8 (similar to that in Fig. 4). As was discussed in section 3.3, overlapping patches quickly annihilate, so the maximum possible number of patches that could form would have to be less than 500 8 ≈ 60. In Fig. 16 we show the average number of patches formed as a function of L scale . We see that for about 3 orders of magnitude, 0.1 < L scale < 10, the number of patches formed is about 18, (around 25% of the maximum possible number).

For large values of L scale fewer patches form, presumably due to the mild gradients causing the dynamics to be similar to the spatially uniform case where essential-extinction results. Finally, for L scale much less than the scale of the dispersal distance (σ 2 = 1) we also see fewer patches form. This can be attributed to the convolution process smoothing out the fine scaled spatial features of g(u 0 (x)), in effect leaving an almost uniform spatial density. Experimenting with increasing the fixed point threshold to 10 -6 or for example running a fixed number of 100 iterations, did not seem to appreciably alter the number of patches formed.

• • • • • • • • • 0.01 0.05 0.1 0.5 1 5 10 50 100 5 10 15 20 L scale avg. number of patches formed

Discussion

In this paper, we studied nonspreading solutions for the integro-difference equation (1) where the growth function g(u) exhibits a strong Allee effect and overcompensation. The nonspreading solutions take forms of stable equilibrium solutions vanishing at ±∞ and solutions oscillating in densities and spatial ranges. Such nonspreading solutions exist in a solid region in parameter space. In a large habitat, patch formation can occur with each patch essentially formed by a nonspreading solution. Our results show that single species model (1) with constant parameters can have very rich nonspreading population dynamics.

Both a strong Allee effect and overcompensation in population growth are necessary to produce biologically meaningful nonspreading solutions in a homogenous environment. Population densities above the Allee threshold generate a forward pushing force. Overcompensation tempers the strength of the pushing force from regions with a high population density. Likewise, a backward pulling force is created from regions where the density is below the Allee threshold.

Nonspreading solutions emerge when there is a balance between the forward pushing and backwards pulling forces in the long run. In the absence of overcompensation, model (1) with a strong Allee effect can have a nonspreading solution if and only if 1 0 [g(u)u]du = 0, so that such a nonspreading solution exists only in a region with measure zero in the parameter space, and thus it is not robust. However, as indicated in the bifurcation diagram Fig. 13, with both a strong Allee effect and overcompensation, there is a solid parameter region (blue) in which nonspreading solutions exist. In this region, 1 0 [g(u)u]du can have any sign. The bifurcation diagram Fig. 8 shows various patterns of nonspreading solutions from steady states, period-two, and several levels of period doubling when a varies while other parameters are fixed.

It should be noted that η, the kurtosis of the dispersal kernel, and initial data also play important roles in developing nonspreading solutions. In Fig. 10, there are no nontrivial nonspreading solutions for η < 2, period-two nontrivial solutions exits on a relatively small interval near η = 2 and for relatively large η there is a stable nonspreading equilibrium. In Fig. 11 and 12, the formation of nonspreading solutions depend on the amplitude and support of initial data. In addition to the growth function used within this paper, the PhD Thesis of Otto [48] demonstrates that nonspreading solutions can form with a variety of other forms of growth functions.

To examine the interaction between nonspreading patch like solutions we considered initial data in the form u 0 (x) = u p x + d 2 + u p x -d 2 where u p (x) is a single patch. In Fig. 15 we were able to show that with sufficient separation these two patch solutions are able to persist for biologically meaningful lengths of time. For example, with d = 8 the population persists for more than 500,000 generations with the parameter values used.

Nonspreading solutions provide a basis for the development of patch formation. For a large habitat, separate patches can emerge from perturbations in a relatively constant population, with each patch basically a nonspreading solution as is demonstrated in Fig. 17. Patch formation is weakly sensitive to the length scale of correlations in the initial distribution as shown in Fig. 16, with patch formation being favored by length scales on the order of the dispersal distance. Correlation lengths much larger or smaller than the dispersal distance result in less patch formation.

We also found that growth parameters giving rise to essential extinction in the non-spatial model, can actually experience population spread and population growth in the spatial model. This is consistent with Vortkamp et al. [57], who using a spatially discrete 2-patch model, demonstrated that essential extinction could be stabilized by an out of phase rescue effect. We will save further investigation of this phenomena for future work.

Our result on nonspreading solutions contrasts with that of Sullivan et al. [55]. Working on (1) with a truncated Ricker's function for population growth, they found that fluctuating spreading speeds can occur as a result of a combination of a strong Allee effect and overcompensation.

The scaled growth function g(u) given by (1) with the carrying capacity 1 has two parameters describing the Allee threshold and strength of overcompensation, respectively. If the carrying capacity of the truncated Ricker's function is scaled to 1, there is no parameter controlling the strength of overcompensation. Therefore the growth function used in this paper is more flexible than the truncated Ricker's function considered in [55]. Observe that in Figure 13, there also exists a solid region (yellow) where spreading solutions exist. Depending on the parameters, model (1) possesses spreading and nonspreading dynamics as well as extinction.

The fact that a single mathematical equation can admit such qualitatively divergent output as spreading solutions, nonspreading solutions, and extinction is intriguing. The possibility of nonspreading solutions is particularly interesting because it suggests a new way to connect the widely employed modeling framework of integro-difference equations to a completely different purpose: the origin and maintenance of ecological boundaries. The factors influencing the location and maintenance of species' spatial distributions, whether patch boundaries on small scales or geographic range boundaries on larger scales, have been the subject of intense interest by ecologists for decades [6,20,54]. The specific biological mechanisms leading to the existence of such boundaries are diverse, but often reflect an interplay between local population dynamics and dispersal. Such dynamics could be related to the oscillating wave fronts observed with this model (see Fig. 4, 5, and [55]). For example, repeated processes of invasion and extinction appear to be important for the maintenance of species' patch boundaries in mixed conifer-hardwood forests [19]. Likewise, Allee effects can contribute to the existence of geographic range boundaries in some insect systems with short dispersal distances [41,50].

Identifying the existence of nonspreading solutions in integro-difference equations opens up several additional lines of inquiry for this modeling framework. One such possibility would involve investigations of how large contiguous populations collapse into small patches, either on evolutionary timescales [45] or in connection with the persistence of relictual populations in conservation biology [10,11]. Likewise, future research could examine nonspreading solutions for integro-difference equations operating on a landscape gradient (e.g., temperature, rainfall) that influences population growth rate. Such studies would provide a vehicle for investigating the interplay between biological and environmental processes that can jointly influence the origin and maintenance of geographic range boundaries [20], including the possibility of patchy population structure at geographic range margins [6,17]. Overall, the existence of nonspreading solutions in integro-difference equations suggests the emergence of a welcome new tool for studying diverse phenomena in spatial ecology.

Acknowledgments

The authors wish to thank Professor Frithjof Lutscher whose detailed review allowed us to greatly improve the rigor and clarity of the manuscript.

Appendix A

Determining an efficient and accurate rule, to code in Matlab, for distinguishing spreading solutions from nonspreading for the model studied here is not a trivial matter. As was shown by Sullivan et al. [55] spreading speeds can fluctuate when Allee and overcompensation are simultaneously present. It is however observed that over a sufficient number of generations the average spread speed will converge to a fixed constant. For spreading solutions we would therefore expect that lim n→∞ spatial extent(u 2n (x)) spatial extent(u n (x)) = 2 .

For the 100 by 100 grid of a and r parameter values scanned in Fig. 13 (a total of 10,000 data points) we see three distinct populations if we look at the ratio of spatial extent(u 500 (x)) spatial extent(u 250 (x)) . Namely, the extinct populations, populations where the ratio is clustered around 1, and populations clustered near 2. For extinct populations, we treat 0/0 as 0. The histogram showing this can be seen in Fig. 18.

The spatial extent of u 500 (x) for the populations with extent ratios near 1 and those with a ratio near 2 differ noticeably. For the population whose extent ratio was between 0.5 and 1.5, we see in Fig. 19 the maximum spatial extent is 4.6. For the population whose extent ratio was greater than 1.5 we see the minimum value of the spatial extent is 6 extending all the way to about 250. This justifies the use of the size extent ratio of 1.5 being used as a threshold to classifying solutions which reach 500 iterations without periodicity being detected. iii. Ethics approval: Not applicable. iv. Consent to participate: Not applicable.