Local interactions affect spread of resource in a consumer-resource system with group defense

Integrodifference equations are a discrete-time spatially explicit model that describes the dispersal of ecological populations through space. This framework is useful to study spread dynamics of organisms and how ecological interactions can affect their spread. When studying interactions such as consumption, dispersal rates might vary with life cycle stage, such as in cases with dispersive juveniles and sessile adults. In the non-dispersive stage, resources may engage in group defense to protect themselves from consumption. These local nondispersive interactions may limit the number of dispersing recruits that are produced and therefore affect how fast populations can spread. We present a spatial consumer-resource system using an integrodifference framework with limited movement of their adult stages and group defense mechanisms in the resource population. We model group defense using a Type IV Holling functional response, which limits the survival of adult resource population and enhances juvenile consumer production. We find that high mortality levels for sessile adults can destabilize resource at carrying capacity. Furthermore, we find that at high resource densities, group defense leads to a slower local growth of resource in newly invaded regions due to intraspecific competition outweighing the effect of consumption on resource growth.


Introduction
Integrodifference equations (IDEs) are a modelling framework that describes a population density in continuous space and discrete time by exploring the growth and dispersal processes separately (Kot and Schaffer 1986).They have been successfully used to study the spread dynamics of annual plants (Andersen 1991), populations in a river system (Lutscher et al. 2010), and populations with moving habitats (Zhou and Kot 2011).This approach has also been expanded to consider population interactions such as consumption (Neubert et al. 1995), parasitism (Cobbold et al. 2005), and competition (Li 2018).
In a variety of organisms such as perennial plants, echinoderms (Black and Moran 1991;Tyler and Young 1998), and colonial insects (Hakala et al. 2019), dispersal happens at some stages in their life history, with other stages being more sessile.The IDE framework can be expanded to consider these dynamics by explicitly adding a non-dispersing stage of a population.Such IDE models find that local interactions of these organisms may limit the number of dispersing recruits that are produced, which may lead to a slower spread rate.For example, Cobbold et al. (2005) and Marculis and Lui (2016) found that parasitism of more sessile stages destabilize the spatial distribution of the entire population and reduce the spread rate.In a model of competition of different green crab genotypes (Kanary et al. 2014), an increased sessile adult survival of the entire population leads to an increased spread rate of the top competitor and a decrease in the spread rate of the lower competitor.However, this dynamic has not been previously included in consumerresource systems, which includes both herbivore-plant and predator-prey interactions.
When considering consumer-resource dynamics in organisms with limited movement, group defense mechanisms may allow resource to become more resistant to consumption.These group defense mechanisms reduce consumption 1 3 intensity as resource density increases (Dubois et al. 2003).This behavior occurs in various taxa where adults have limited movement.For example, bees produce social waves that repel predators (Kastberger et al. 2008), and kelp provide habitat for predators of their grazers, which induces cryptic behavior in grazers with subsistence off of kelp detritus rather than active grazing (Karatayev et al. 2021).A previous model without stage-dependent dispersal considered the spatial dynamics of a resource using group defense (Venturino and Petrovskii (2013)), where they found oscillatory spatial distributions at high initial resource densities caused by group defense.The potential for group defense to qualitatively affect dynamical outcomes of interacting species raises the question of how group defense in a sessile stage might affect overall spread given dispersive juveniles.
In this paper, we present and analyze an IDE model of the spatiotemporal dynamics of a consumer-resource system where adults have limited movement and resource present group defense.In "Model", we introduce the model, which is based on the ideas presented by Kanary et al. (2014), and provide a nondimensional version which we will analyze.In "Results", we analyze two features of the spatiotemporal dynamics: the dispersal-induced instabilities of the resource-only system and the spread rate of resource.Finally, in "Discussion", we discuss how these results lead to a further understanding of how local interactions affect the spread of organisms.

Model
In this section, we extend an integrodifference model to consider motile, dispersing juveniles and the local interactions between sessile adult stages.A similar extension was previously considered in Kanary et al. (2014), and a formal construction of this model is analogous to that in Arroyo-Esquivel et al. (2022).Consider a region in 1-dimensional space denoted by Ω ⊆ ℝ .At each time step m and point in space x ∈ Ω , our model follows population densities of con- sumer P m (x) and resource N m (x) populations at reproductive age (hereafter adults).
For each population i = P, N , at each time step m, a fraction i of the adult population survives to the next time step in the absence of consumption.Consumers consume adult resource following a unimodal, Type IV Holling functional response (Andrews 1968).This functional form, denoted by G(P, N), has an attack intensity of the consumer N and a group defense intensity parameter N .With these parameters, the Type IV Holling functional response takes the form: A simple calculation shows that this functional form reaches its maximum consumption intensity at a resource density (1) of 1∕ √ N .At higher resource densities, the consumption intensity decreases.
Juveniles of both populations disperse following a kernel k i , for i = P, N , where k i (x, y) represents the probability density function of a dispersing individual starting in point y ∈ Ω to arrive at point x.Following dispersal, a fraction of those juveniles survive and become adults at the next time step.Consumers produce juveniles that survive to become adults proportional to their consumption intensity with a factor P .Resources produce juveniles by a constant percapita factor R, where R > 1 −  N (i.e., more juveniles are produced than adults die) for population persistence.The fraction of newly setting resource juveniles that survive to become adults depends on local consumer and resource densities.Consumers consume settling resources before they can provide with group defense mechanisms with a constant intensity S .Local resources further limit resource settlement through intraspecific competition with intensity .This makes the resource density have a carrying capacity proportional to 1∕ .
These assumptions lead to the system of equations: A summary of the parameter's meanings and dimensions is found in Table 1.To simplify our analysis, we first nondimensionalize the model.We use the same nondimensionalization as in Arroyo-Esquivel et al. (2022).For each m, let p m = S P m , n m = N m .Then, if p = P ∕ , and n = N ∕ S , and = N ∕ 2 , our nondimensional version of the model is (2) Note that we have also changed the indices of i and k i in order to preserve clarity.
In the case the dispersal dynamics in Eq. 3 are disregarded (by choosing k i (x, y) = (x − y) the Dirac delta), our system reduces to that studied in Arroyo-Esquivel et al. (2022): The fixed points of Model 4 and their stability are relevant for the analyses we will perform in "Results".Here, we present those results, also summarized in Table 2. Model 4 has four fixed points: an unstable extinction fixed point (0, 0), a resource-only fixed point (0, n * ) where n * is and, when ≠ 0 , two coexistence fixed points (p ∧ , n ∧ ) and (p ∨ , n ∨ ) .The lower coexistence fixed point (p ∨ , n ∨ ) is always unstable, whereas the higher coexistence point (p ∧ , n ∧ ) exchanges stability with the resource-only point (0, n * ) via a transcritical bifurcation at the bifurcation value for p : The resource-only fixed point is stable when consumer conversion intensity is under a given threshold * p (i.e. p <  * p ) and unstable otherwise.In addition, the positive coexistence fixed point becomes biologically infeasible as it becomes stable as p ∧ becomes negative.We can thus (3) (4) say that the positive coexistence fixed point is unstable whenever it is biologically feasible.
In the case the resource-only equilibrium is stable (when  p <  * p ), this stability is global, i.e., all trajectories converge to the equilibrium.In the case there are no stable fixed points (when  p >  * p ), the system converges globally to a quasiperiodic consumer-resource cycle.

Dispersal-induced instabilities
In this section, we explore how dispersal affects the stability of the consumer-resource dynamics by analyzing Model 3. Dispersal can induce instabilities in stable population densities, a mechanism first observed by Turing (1990) and further analyzed by Levin (1974).In the case of IDEs, this mechanism can be analyzed following the linearization process of Neubert et al. (1995).Throughout this section, we assume that the dispersal kernels only depend on the distance between the two points x, y ∈ Ω , i.e., k i (x, y) = k i (x − y) .We can write Model 3 in the same form as the model presented in Kanary et al. (2014): In addition, we will assume that dispersing juveniles do not die or escape the habitat during the dispersal process, i.e., for both i = p, n .Using this assumption, we linearize the sys- tem near a stable equilibrium as follows.Let (p, n) be a stable equilibrium of System 4, then if where ( m , m ) is a small perturbation around (p, n) , we lin- earize the first equation of System 7 as (7) Multiplying the terms around the integral and disregarding higher-order terms yields and a similar equation for m+1 .Then, given J(F) as the Jacobian matrix of a given function F evaluated at (p, n) , our linearized system, in matrix form, is where K(x, y) = diag(k p (x − y), k n (x − y)) and the integral represents rowwise integration.To study how dispersal leads to instabilities in the system, let  p <  * p , where * p is defined by Eq. 6.This makes the resource-only fixed point (0, n * ) to be stable.Thus, let (p, n) = (0, n * ) .Then, the linearized system (Eq.8) is We then take the Fourier transform of Eq. 9.By doing this, our system simplifies to where f corresponds to the Fourier transform of a given function f, i.e., To properly apply the Fourier transform to the functions (x) and (x) , we assume that (x) = (x) = 0 for x ∉ Ω .The matrices A,K,J in Eq. 10 satisfy Decay of ξm () and ηm () for all guarantees the decay of m (x) and m (x) , which would imply stability of the carry- ing capacity equilibrium.The matrix A+KJ has a triangular form, and thus, the eigenvalues are If both populations disperse following a Laplace kernel with mean dispersal distance 1∕a i for i = p, n, then the Fourier transform of these kernels is Note that 0 ≤ ki () ≤ 1 for all , which implies that if  p <  * p , then 0 <  1 < 1 and  2 < 1 for all .For 2 , for any R,  n > 0 , the inequality  2 < −1 does not have a real solu- tion.This implies that dispersal of a Laplace kernel will not induce instabilities in a resource-only state.
If we choose instead a double-gamma distribution, for i = p, n , then their Fourier transform is Similar to in the case of Eqs. 15 and 17 satisfies ki () ≤ 1 , which implies  1 < 1 and  2 < 1 for all .Equation 17 has a global minimum of − 1 8 .If kp () = −1 8 , then we find that  1 < −1 is satisfied when , and thus, instabilities will only be caused by dispersal for high local resource mortalities.We can compare that, at low local mortalities of resource, the eigenvalue 2 is almost unchanged as the frequency changes (Fig. 1a), whereas at high local mortalities, the eigenvalue has a wider range of change and crosses ( 13) the threshold of = −1 (Fig. 1b).This suggests that disper- sal of resource does not affect its stability when most of the reproductive adults can survive more than one time step.
The spatial pattern formation presented by these dynamics is in Fig. 2.Even in the presence of instabilities caused by the dispersal of resources, consumers are not able to invade.This shows that although resource density is varying, a low consumer conversion rate (  p <  * p ) makes it impossible for consumers to invade and have any influence over the resource population, making this system essentially a resource-only system.

Spread rate of resource
In a general integrodifference framework, the spread rate of a population is calculated by analyzing when the extinction equilibrium of the system (0, 0) becomes unstable (Zhou and Kot (2013)).This also works in the case of a single-population dynamics system with sessile stages (Cobbold and Stana (2020)).In the case of System 3, the extinction equilibrium of resource (0, 0) is always unstable, which implies that the resource is always able to invade when rare.Instead of explicitly calculating the Fig. 1 Values of 2 in Eq. 13 as the Fourier transform frequency varies, when k n follows the double-gamma distribution kernel (Eq.16).In these figures, we use R = 20 , and a n = 1 with two values for resource adult survival: a n = 0.9 and b n = 0.01 .Note that, for high resource mortalities (subfigure b), some frequencies of the Fourier transform produce an eigenvalue  2 < −1 , which induces instabilities caused by dispersal spread rate, we numerically estimate the time it takes for the population to reach a specific population density at a given point in space.
To do this, let Ω = [−L∕2, L∕2] for habitat length L. In addition, let the initial conditions be a constant consumer density p 0 (i.e., p 0 (x) = p 0 ) and the resource at carrying capacity at a single point at 20% of the length of the habitat (i.e., n 0 (x) = n * (x − a) , where (x) is the Kronecker delta, and a is the point that represents 20% of the length of Ω ).We then calculate the time it takes for the resource to reach a population density of 80% its carrying capacity at 80% of the habitat length, i.e., we find the time M that satisfies where b is the point that represents the 80% of the length of the habitat.An example of this procedure is in Fig. 3.In this case, M = 28 .We then explore how changing different parameters of the model makes these transient times vary.
The results of these numerical experiments are in Fig. 4. Intuitively, as more consumers are present in the environment (higher initial consumer density p 0 and consumer sur- vival p ), the time to spread increases, and as more resource juveniles are produced (higher R (h)) and disperse further (higher a n (i)), the time to spread decreases.Notice however Fig. 2 Distributions of consumers p t (x) and resources n t (x) after 1000 time steps with an initial distribution being a random perturbation of the uniform distributions p 0 (x) = 0 and n 0 (x) = n * .In these simulations, we use a doublegamma kernel (Eq.16) with a n = 1 and a p = 5 .The other parameters of the model are p = 0.8, N = 0.1, = 1, R = 20, p = 0.7 , and * p with two values for consumer adult survival: a n = 0.9 and b n = 0.01 .In the case of low consumer mortality (a), the system is unable to escape the fixed point and converges to the uniform distribution of the fixed point (0, n * ) .In the case of high consumer mortality (b), the dispersal-induced instabilities cause the resource to fluctuate through space.Even in this case, the low conversion efficiency of p consumers makes the consumer population unable to invade that too high of a resource dispersal a n leads to an increase in the time it takes for the population to reach a high density.This is caused due to a loss of dispersing individuals at the edges of the domain.More surprisingly, increasing adult resource survival (higher n ) and increasing the attack intensity ( n ) lead to an increase in this time, and consumer conversion intensity ( p ) and group defense intensity ( ) have a minimal impact.
To explain these results, note that a higher adult resource survival n leads to a higher value of n * in Eq. 5.This higher value of n * takes longer to be reached, thus making the conditions of Eq. 18 take longer to be satisfied.Although group defense intensity itself did not have an impact over the time it takes for the population to reach carrying capacity on the other of the domain, group defense explains why the consumption might not have a big impact over the time to spread, which is the case when there is no group defense at offspring production intensities orders of magnitude smaller than when there is group defense (Fig. 5).This is because group defense makes consumption less important at high resource densities, which means the source population at the left of the habitat is not impacted by consumption and allows to source new juveniles that will eventually Fig. 3 Distributions of consumers p t (x) and resources n t (x) at a initial setup and b after 30 time steps.In these simulations, we use a Laplace kernel (Eq.14) with a n = 1 and a p = 5 .The other parameters of the model are L = 10, p = 0.8, N = 0.1, = 1, R = 20, p = 0.7 * p , and n = 0.9 Fig. 4 Number of time steps M required for the population to reach 80% of its carrying capacity at 80% of the habitat Ω (Eq.18) when varying a single parameter in Model 3. Unless it's the parameter being changed, in these simulations, we use a Laplace kernel (Eq.14) with a n = 1 and a p = 5 .The other parameters of the model are L = 10, p = 0.8, N = 0.1, = 1, R = 20, p = 0.1 , and n = 0.9 overcome consumption at a similar density.A similar argument explains the low impact of consumer conversion intensity p over the spread time.
To make the value M independent of the parameters of the model, we repeat the analysis with an arbitrary threshold that does not depend on any parameters.By choosing such threshold to be equal to 1, in Fig. 6, we observe that resource survival ( n ) does not have a big impact on spread speed, while most of the other parameters have a similar behavior to the one observed in Fig. 4. The most likely unintuitive difference is that consumer survival ( p ) does not affect in any way the ability of the resource to spread to the other side of the domain.However, this makes sense as consumption at low densities is weak, and thus, the main factor affecting spread is how much of the resource is able to survive once settled.This also explains the small decrease in the time M as resource survivability becomes high enough.

Discussion
In this paper, we explored how local interactions of sessile organisms in a consumer-resource system affect the spread rate of resource.Two main results arise from this exploration: sessile resource adults stabilize the spatial distribution of resource, and group defense leads to a slower spread rate.
We find that, when the consumers cannot invade the resource, and when most adults survive to the next reproduction period (high n ), these sessile adults stabilize the distribution of resource and prevent the resource carrying capacity to be destabilized by dispersal.This destabilization required a fat-tailed kernel, which leads to accelerated invasions (Kot et al. 1996).These stability results also provide more evidence to the argument that increased dispersal leads to a negative correlation between spatial stability and synchrony in population densities between patches (Abbott 2011).
In the case of resource spread when rare, we found no effect of increasing the group defense intensity ( ) in the spread rate of the resource.However, group defense allows.There is evidence that would suggest group defense should have a negative impact on the spread rate.For example, the decrease in consumption could lead to higher intraspecific competition, which affects juvenile survival (Morin 1986;Wilson 1989) and spread (Oricchio and Dias 2020).However, none of these studies looked at this question in the context of organisms presenting group defense.Another way that group defense might decrease spread rate, not modelled here but biologically feasible, is if resources' energy investment in group defense reduces energy availability for reproduction (Sasmal and Takeuchi 2020).Our model suggests these impacts are compensated by the reduced effect of increased consumption, which does have an effect when group defense is not included in the model ( = 0 ), which explains why has no effect on the spread of the population.
Another feature of our model is the implicit inclusion of stage structure in the spatial dynamics.The importance of stage structure in dispersal dynamics was first observed by Hastings (1992) using spatially explicit models with continuous space and continuous time and with discrete space and discrete time.In the case of continuous space and discrete time, previous analyses that expand the IDE framework to consider stages with limited movement also implicitly included stage Fig. 5 Number of time steps M required for the population to reach 80% of its carrying capacity at 80% of the habitat Ω (Eq.18) when varying the consumer offspring production in Model 3 in the case where there is no group defense.In these simulations, we use a Laplace kernel (Eq.14) with a n = 1 and a p = 5 .The other parameters of the model are L = 10, p = 0.8, N = 0.1, = 0, R = 20 , and n = 0.9 Fig. 6 Number of time steps M required for the population to reach a density of 1 at 80% of the habitat Ω (Eq.18) when varying a single parameter in Model 3. Unless it's the parameter being changed, in these simulations, we use a Laplace kernel (Eq.14) with a n = 1 and a p = 5 .The other parameters of the model are L = 10, p = 0.8, N = 0.1, = 1, R = 20, p = 0.1 , and n = 0.9 structure (Kanary et al. 2014;Veit and Lewis 1996).We find a similar general result to those previous analyses, where local interactions in the stages with limited movement affect the spread rate of the population.In our model, intraspecific competition slows down spread, in contrast to competitive systems, where a high survival of adults promotes the spread of the top competitor (Kanary et al. 2014), and analogous to single population dynamics where high mortality of stages with limited movement can lead to an Allee effect which slows spread rates (Veit and Lewis 1996).
As with any model, we made a number of simplifying assumptions in our model.First, we only consider the case where consumers resource upon adult resource or settling juveniles.Growth of species with limited movement as adults such as urchins (Allen 2008) and tunicates (Olson and McPherson 1987) has been linked to the consumption of their dispersing larvae.We suspect dispersing larvae consumption will reduce the impact of local interactions and give more emphasis on the dispersal dynamics.Second, we assume that the dispersing individuals are the juveniles.This assumption does not capture populations where larvae have limited movement ability compared to their adult stages such as a consumer-resource interaction between dragonflies and frogs (Caldwell et al. 1980).We speculate the model that captures those dynamics would have a similar structure to this one, which would imply that juvenile interactions would be the main constraint on of spread dynamics.
Another limitation is that our environment is spatially homogeneous.In reality, spatial heterogeneity may lead to different dynamics than the ones observed in our model.In our model, the only factor that produces habitat heterogeneity for resource is the distribution of consumer population.However, other potential factors of heterogeneity not accounted by our analysis are substratum topography (Erlandsson et al. 2005;Köhler et al. 1999) and resource availability (Grabowska and Kukliński 2016).These factors can be included in our model with spatially variable survival or reproductive functions.This could render our problem intractable, which would require numerical analysis to be well explored.
Finally, we assumed both consumer and resource have limited movement as adults.However, by setting the proportion of sessile adults that survive ( i for i = p, n ) equal to 0, our model allows only one of the two species to have limited movement as adults.As seen in the dispersal-induced instabilities ("Spread rate of resource"), this is a sufficient condition for instabilities of resource at carrying capacity.
In conclusion, in a consumer-resource system, local interactions between sessile adults are key to determining the ability of resource to spread by limiting their production of offspring through consumption.Similar results were obtained when modelling invasive algae, where the consumption of the substrate in the soil slowed the spread rate of the algae (Britton-Simmons and Abbott 2008).These observations contrast with those seen in competitive models, where competition acts on a more regional scale by allowing the coexistence of competitors in space (Allen et al. 1996) or stopping the invasion front of the higher competitor (Kanary et al. 2014).These models exemplify the use of the IDE framework in a wider range of interactions between species such as perennial plants and animals with limited movement.

Table 2
Summary of the stability results of the fixed points of Model 4 p (p ∧ , n ∧ ) Unstable for  p <  * p (Eq. 6) and stable (but biologically infeasible) for  p >  * p (p ∨ , n ∨ ) Always unstable