Nonspreading solutions and patch formation in an integro-difference model with a strong Allee effect and overcompensation

Previous work involving integro-difference equations of a single species in a homogenous environment has emphasized spreading behavior in unbounded habitats. We show that, under suitable conditions, a simple scalar integro-difference equation incorporating a strong Allee effect and overcompensation can produce solutions where the population persists in an essentially bounded domain without spread despite the homogeneity of the environment. These solutions are robust in that they occupy a region of full measure in parameter space. We develop orbit diagrams showing various patterns of nonspreading solutions from stable equilibria, period-two, to higher periodicity. We show that from a relatively uniform initial density with small stochastic perturbations, a population consisting of multiple isolated patches can emerge. In ecological terms, this work suggests a novel endogenous mechanism for the creation of patch boundaries.


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, integrodifference models have been used to predict changes in gene frequency (Lui 1982a(Lui , b, 1983;;Slatkin 1973;Weinberger 1978), and characterize species' spatial dynamics (Hardin et al. 1988a(Hardin et al. , b, 1990;;Hastings and Higgins 1994;Kot and Schaffer 1986;Kot 1989Kot , 1992Kot , 2001;;Neubert et al. 1995).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 (Kot et al. 1996) and cases where one or more species can (or cannot) outrun the pace of environmental change (Hu et al. 2020;Li and Wu 2020;Zhou and Kot 2011).The reader is referred to the monograph by Lutscher (2019) 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. (2001) 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 (Fagan et al. 2009).Related results for gap crossing in integro-difference equations can be found in Musgrave et al. (2015).
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 (Allee et al. 1949;Ashman et al. 2004;Burd 1994;Cadre et al. 2008;Calabrese and Fagan 2004;Courchamp et al. 2008Courchamp et al. , 1999;;Davis et al. 2004;Dennis 1989;Groom 1998;Larson and Barrett 2000;Moeller 2004;Parker 2004;Stephens and Sutherland 1999), and they have been studied in connection with integro-difference equations in the context of spatial spread (Lewis and Kareiva 1993;Wang et al. 2002).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 (Calabrese and Fagan 2004;Groom 1998), 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 (Lynch et al. 2014;Rhainds and Fagan 2010).
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. (1986) 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., Barraquand et al. 2014;Framstad et al. 1997).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(x − y) .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 overcompensa- tion, 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 1982;Li et al. 2009).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. 2018(Bourgeois et al. , 2020;;Li et al. 2009).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 overcom- pensation).In this case, the spreading speed is the unique speed of traveling waves connecting zero and the carrying capacity (Lui 1983), and the sign of the wave speed is the same as that of ∫ 1 0 [g(u) − u]du (Wang et al. 2002).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. 2017).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 Sullivan et al. (2017), 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 (1) 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. (2020).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 sensi- tive 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.
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 (Schreiber

(a) (b)
Fig. 1 The sensitivity of the maximum value of g a,r (u) to the parameters a and r

2003
).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 Wang et al. 2002).While this is not neces- sarily 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 (Kot et al. 1996).To model dispersal with varying kurtosis, we will use the generalized Gaussian distribution (Nadarajah 2005) 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 expo- nentially 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 hav- ing 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 cal- culations.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/ glott o01/ theor etical-ecolo gy.git.Numerical experimentation showed us that decreasing or the clipping threshold did not alter results.For example, Fig. 8 is recreated using a clipping threshold of 10 −5 and = 0.0025 with identical results.
Fig. 2 The graph of g a,r (u) and y = u showing the effects of varying a and r

Nonspreading solutions
In contrast to integro-difference equations with Allee or overcompensation effects considered separately (Li et al. 2009;Lui 1983;Wang et al. 2002), 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.
In Figs.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 Fig. 4 Solution with parameters a = 0.62, r = 8, = 5, p 0 = 0.9, w 0 = 4 .In part (a) the blue curve is u 0 (x) , and the gray curves are the transients u 1 through u 15 and the black is u 16 through u 100 .In (b) we see the spatial extent converging to that of the stable equilibrium sample points are geometrically spaced).As can be seen in Figs.4-6, 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 behavior 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 (Figs. 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.

Two-parameter bifurcations
The Matlab codes used in this section can be found in the folder 2parameter_bifurc on https:// github.com/ glott o01/ theor etical-ecolo gy.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 param- eter 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 a − r 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) (b)
Fig. 5 Solution with parameters a = 0.585, r = 8, = 5, p 0 = 0.9, w 0 = 4 .In part (a) the blue curve is u 0 (x) , the gray curves are the transients u 1 through u 16 , the red are the odd indexed iterations u 17 through u 99 , and the black are the even indexed iterations u 18 through u 100 .In (b) we see the spatial extent oscillating with period of length two 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.
Next we examine the dependence on initial conditions by varying the parameters w 0 and p 0 in model (1) for fixed a and r.Here we present the case where a = .61,r = 8, = 5 but the results are qualitatively representative of what is seen for other parameter choices that have nonspreading solutions based on our extensive simulations.Three phenomena are seen to occur: extinction, the formation of a single unimodal population "patch", or the formation of two spatially disjoint unimodal "patches" of population.If a unimodal equilibrium occurs, its shape is independent of initial conditions and it appears similar to that in Fig. 4.
The phenomena of two patches forming are shown in Fig. 7.We see initially the density near x = 0 experiences high growth, overcompensation causes this to subsequently fall below the Allee threshold, thus effectively separating the left and right sides of the population.In Fig. 14 we show a bifurcation diagram for p 0 and w 0 .Blue is used for a single unimodal patch, gray for extinction, and red for the formation of two unimodal patches.Predictably extinction occurs if the support of the initial domain is too small (small w 0 ), or if the initial density is too small (small p 0 ).For values of p 0 which correspond to high growth (roughly 0.8 ≤ p 0 ≤ 1.2 ) and a length of support comparable to the spatial extent of the unimodal equilibrium (roughly 3 ≤ w 0 ≤ 7 ) we see a single patch emerge.For larger values of w 0 we see two-patch solutions emerge.
To create Fig. 14 we iterated until a fixed point condition was met, namely that the maximum absolute difference between subsequent density curves was less than 10 −5 .Once the fixed point condition was reached, the type of equilibrium was determined by integrating the population density.We found that the total population was very nearly an integer multiple of the total population of the unimodal equilibrium which is 3.695.Only multiples of 0, 1, or 2 were observed.

Interaction of two patches
The Matlab codes used in this section can be found in the folder two_patch on https:// github.com/ glott o01/ theor eticalecolo gy.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 (Schreiber 2003), 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 form u 0 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 super-exponentially (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.

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/ glott o01/ theor eticalecolo gy.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 [0.8, 1.2].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 the "Interaction of two patches" section, 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 num- ber 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 pos- sible 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 leav- ing 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.

Discussion
In this paper, we studied nonspreading solutions for the integro-difference Eq. ( 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 den- sities 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 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 exist on a relatively small interval near = 2 and for relatively large there is a stable nonspreading equilibrium.In Figs.11 and 12, the formation of nonspreading solutions depends on the amplitude and support of initial data.In addition to the growth function used within this paper, the PhD Thesis of Otto (2017) 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 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. (2020), 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 phenomenon for future work.
Our result on nonspreading solutions contrasts with that of Sullivan et al. (2017).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 Sullivan et al. (2017).Observe that in Fig. 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 (Brown et al. 1996;Gaston 2009;Strayer et al. 2003).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 Figs. 4, 5, and Sullivan et al. 2017).For example, repeated processes of invasion and extinction appear to be important for the maintenance of species' patch boundaries in mixed conifer-hardwood forests (Frelich et al. 1993).Likewise, Allee effects can contribute to the existence of geographic range boundaries in some insect systems with short dispersal distances (Lynch et al. 2014;Rhainds and Fagan 2010).
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 (Nekola 1999) or in connection with the persistence of relictual populations in conservation biology (Channell and Lomolino 2000a, b).Likewise, future research could examine nonspreading solutions for integrodifference 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 (Gaston 2009), including the possibility of patchy population structure at geographic range margins (Brown et al. 1996;Fortin et al. 2005).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.

Appendix
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. (2017) 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 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.

Fig. 3
Fig. 3 Generalized Gaussian distributions with standard deviation equal to one for various values of

Fig. 14
Fig. 14 Initial condition parameters resulting in a single patch are shown in blue.Red represents parameters giving rise to two patches.Gray represents parameters resulting in extinction.Values used for other parameters are a = .61,r = 8, = 5

Fig. 16
Fig. 16 The mean number of patches formed as a function of the correlation length scale of the stochastic initial condition.The black bars are the 90% confidence intervals.Note the x-scale is logarithmic

Fig. 17
Fig. 17 Typical plots of initial conditions (blue) and resulting patch formation (red) at short, intermediate, and long length scales

Fig. 18 Fig. 19
Fig.18Histogram of the ratio of the spatial extent for u 500 to that of u 250 for the parameters scanned in Fig.13.Extinct populations are binned in x = 0