Characterizing Long Transients in Consumer–Resource Systems With Group Defense and Discrete Reproductive Pulses

During recent years, the study of long transients has been expanded in ecological theory to account for shifts in long-term behavior of ecological systems. These long transients may lead to regime shifts between alternative states that resemble the dynamics of alternative stable states for a prolonged period of time. One dynamic that potentially leads to long transients is the group defense of a resource in a consumer–resource interaction. Furthermore, time lags in the population caused by discrete reproductive pulses have the potential to produce long transients, either independently or in conjunction to the transients caused by the group defense. In this work, we analyze the potential for long transients in a model for a consumer–resource system in which the resource exhibits group defense and reproduces in discrete reproductive pulses. This system exhibits crawl-by transients near the extinction and carrying capacity states of resource, and a transcritical bifurcation, under which a ghost limit cycle appears. We estimate the transient time of our system from these transients using perturbation theory. This work advances an understanding of how systems shift between alternate states and their duration of staying in a given regime and what ecological dynamics may lead to long transients.


Introduction
One of the goals of mathematical modeling of ecological systems is to understand the fate or long-term dynamics of such system. The main method to study such fate has been through the analysis of the attractors in a model (Ives and Carpenter 2007). Recent years have seen an increase in the interest of understanding non-attractor dynamics (hereafter transients) of the models, especially those that resemble an attractor for a long period of time (hereafter long transients) (Hastings et al. 2018). Long transients have gained recognition as a theoretical tool to better describe population dynamics by allowing the study of dynamics that occur in a more biologically relevant timeframe (Morozov et al. 2020). In addition, an understanding of long transients can inform conservation and natural resource management goals. For example, identifying that a positively valued long-term behavior observed in nature is actually a long transient and what causes it can guide management to prolong it (Francis et al. 2021).
Long transients often appear in the presence of a "small" (close to zero) parameter in the model (Morozov et al. 2020). One of the main challenges of identifying long transients is identifying such a small parameter, which may be a function of the biologically reasonable parameters, and thus may not be easily interpretable. For example, in ghost attractors, this small parameter is the difference between a bifurcation parameter and its bifurcation value (Morozov et al. 2020). While varying the parameter past such bifurcation leads to the destruction of an attractor, small differences the transient dynamics will resemble the attractor. In crawl-by attractors, the small parameter is determined by the degree to which the trajectory of the system is parallel to the stable manifold of a saddle node equilibrium at a given time (Morozov et al. 2020). In this case the system will behave similarly to such a stable manifold for a prolonged period of time before the unstable part of the trajectory leads to a change in the system behavior.
One behavior that has been demonstrated to lead to long transients in consumerresource systems is group defense (Venturino and Petrovskii 2013). Group defense is a behavior where a resource population reduces the risk of individuals being predated by protecting each other. This behavior occurs in diverse animal taxa, which produce early-warning signals to detect predators, as is the case of colonial spiders (Uetz et al. 2002), birds, (Robinson 1985), and mammals (Ebensperger and Wallem 2002). This behavior also occurs in producers such as kelp, where high densities of kelp lead to an increase in predators of kelp grazers, which induces cryptic behavior on such grazers and thus reduces grazing intensity (Karatayev et al. 2021).
Group defense transients might also depend on lags in population growth caused by discrete reproductive pulses. In some taxa that exhibit group defense, adult stages of the population may reproduce in discrete, seasonal pulses, such as is the case of kelp (Karatayev et al. 2021) or bees (Kastberger et al. 2008). This can provide individuals to a population decades after stressful events which cause population declines, such as competitive exclusion of pioneer species in tropical rain forests (Dalling and Brown 2009), or extreme weather events in phytoplankton (Ellegaard and Ribeiro 2018).
In this paper we characterize the long transients in a consumer-resource with both group defense and reproductive pulses. We first construct the model that describes a consumer-resource interaction where the resource exhibits group defense and has discrete reproductive pulses. Then, to illustrate the long transients present in this model, we identify a small parameter that describes each of the transients (crawl-by and ghost attractor), and we use this parameter to calculate the time the system remains in this long transient (hereafter transient time). Finding approximations for these parameters and transient times provides biological insight into how these long transients may arise in natural systems with the modeled dynamics. We conclude this paper with a discussion of these results and their biological implications.

Model
In this section we construct a consumer-resource model with group defense and discrete reproductive pulses. We previously explored a spatial, non-smooth version of this model to understand spread of kelp being grazed by urchins . We consider the dynamics of adult consumer P and adult resource N densities through time. Adults of population i = P, N experience a natural mortality at a rate d i . In addition, consider that consumers consume resource following a unimodal Type IV Holling functional response that represents group defense with a decline in consumption at high resource densities (Andrews 1968). We let γ N be the attack rate of the consumer, and the maximum per-capita resource consumption occurs when N = 1 √ σ N . Reproduction and recruitment of juvenile stages occur at discrete points in time. We model this recruitment as an impulsive differential equation. Let t = m be the periods at which the offspring recruit to the population. The number of consumer recruits is proportional to the amount of resource consumed at time t = m with proportionality constant γ P . Resource produce a per-capita number R of recruits. We assume that R > 1 − exp(−d N ) in order to have a self-replenishing resource in the absence of consumers. For predation, a fraction of those offspring survive consumption with a probability following an exponential distribution with mean 1 γ S . Resource offspring also survive intracompetition from adults with carrying capacity proportional to 1 β . Then, given P − m+1 as the density of consumers before the pulse and P + m+1 its density after the pulse (with analogous notation for resource, N − m+1 and N + m+1 ), the dynamics of the adult consumer and resource populations satisfy the following system of impulsive differential equations: (1) We next transform Model 1 into a discrete-time model. We can rewrite the continuous part of the Model 1 as (2) Following the derivation of Cui et al. (2016), we discretize System 2 as By taking δ P = exp(−d P ), δ N = exp(−d N ) and using P − m+1 = P m and N − m+1 = N m as described in System 1, we arrive the following discrete-time model: To simplify our analysis, we will study a nondimensional version of the model. For each m, let p m = γ S P m , n m = β N m . Then, if γ p = γ P /β, γ n = γ N δ P /γ S , σ = σ N /β 2 , our nondimensional version of the model is Note that we have also changed the indices of δ i and k i in order to preserve clarity.

Analysis and Results
In this section we characterize the dynamics of Model 5 and its potential for long transient dynamics. We identify two different classes of long transients, a crawl-by transient around the extinction of resource and another around the carrying capacity of the resource, and a ghost consumer-resource cycle. To illustrate the transients identified and test the accuracy of our analytical approximations, we also characterize all transients numerically by iterating the logarithm of Model 5 in Julia, where the used, fixed parameters and initial conditions are specified as relevant to each analysis below. Based on preliminary numerical simulations in double-precision floating-point numbers, we find that iterating the logarithm of Model 5 instead of the original model prevents numerical instabilities potentially occurring at long-transients near zero by reducing the range of the derivative near zero. The source code for these simulations can be found in https://github.com/jarroyoe/characterizing-transients.
We analytically derive approximations for the transient time of the crawl-by transients using perturbation theory, while we numerically analyze the ghost attractor transient time by regressing the transient time using a power law, which we describe further in Section 3.b. Before we characterize these long transients, we first analyze the equilibria of the model. This model has up to four biologically relevant fixed points: a resourceonly carrying capacity equilibrium (0, n * ), an unstable extinction equilibrium (0, 0), and two possible unstable coexistence saddle equilibria ( p ∨∧ , n ∨∧ ). We assume that the carrying capacity of resource is greater than the density at which consumption growth is its highest, i.e., n * > 1/ √ σ , such that group defense is relevant to resource populations below carrying capacity. Under this condition, the equilibria (0, n * ) and ( p ∧ , n ∧ ) go through a transcritical bifurcation at In this case, the equilibria (0, n * ) is stable for γ p < γ * p and unstable for γ p > γ * p , and the equilibrium ( p ∧ , n ∧ ) is unstable for γ p < γ * p and is not in the first quadrant (i.e., R 2 + ) for γ p > γ * p . See Appendix A for the expressions of these equilibria and their stability. This analysis allows us to better understand the nature of the transients we have identified.

Crawl-by Transients
Although the coextinction equilibrium is a saddle in the n-direction (which implies that n will stay above 0), System 5 can resemble a system where the resource is extinct for a long period of time when consumer density is high (Fig. 1). This is an example of a long crawl-by transient. We determine how prevalent this behavior is in the following theorem, proven in Appendix B.
Theorem 1 Let ε 1. If p 0 is of order ε −1 and n 0 of order 1, then System 5 goes through a crawl-by transient at the extinction of resource n = 0. Recovery of resource will begin after a time of approximately In Fig. 2 we show that Eq. 7 is a reasonably close estimator of the observed point of recovery, especially for large values of the initial consumer density p 0 , where the consumer density is an order of magnitude higher than the initial resource density n 0 . We will show in the following theorem that the long-term dynamics shown in Fig. 1 did not depend on the initial conditions of the model.

Theorem 2 System 5 has a compact, connected global attractor in the first quadrant
See Appendix C for the proof of this theorem. Theorem 2 implies that, when γ p < γ * p , System 5 will go towards carrying capacity of resource and extinction of consumers. On the other hand, when γ p > γ * p , there are no stable fixed points in the first quadrant. Thus, Theorem 2 implies the existence of a nonlinear attractor, which we can describe based on numerical observations, as shown in Fig. 3.
When γ p > γ * p , the resource population is able to reach a maximum density of carrying capacity and stay there for a prolonged period of time (Fig. 3). However, after the consumer reaches a high enough density, the resource population collapses and passes through a transient extinction phase. This cycle repeats itself through time, but at each repetition, the amplitude of consumer density varies. We hypothetize that this variation in amplitude is caused by the system having a long periodicity. In addition, increasing γ p increases the period between each oscillation. This is consistent with the implication from Theorem 1 that a higher consumer density causes the resource to stay around the extinction equilibrium for a longer period of time.  Figure 3 also shows that the system can stay around the resource-only equilibrium for a prolonged time. We approximate this time in the following theorem, proven in Appendix D.
Theorem 3 Let γ p > γ * p , where γ * p is defined by Eq 6 and let Then, if ( p 0 , n 0 ) = (ε, n * − ε) for 0 < ε 1, System 5 goes a crawl-by transient at the resource-only equilibrium n = n * . resource will start decaying after a time of approximately In Fig. 4 we show that this expression is a reasonably close approximation of the time it takes for the consumer to escape extinction across a wide range of orders of magnitude for the initial consumer density.

Ghost Attractors
Theorem 2 ensures that, when γ p < γ * p , System 5 will converge to the stable equilibrium (0, n * ). However, when γ * p − γ p 1, this convergence can take a significantly longer time, as can be seen in Fig. 5. Before the system reaches the equilibrium, the dynamics resemble pseudo-oscillations similar to those observed in Fig. 3 when γ p > γ * p . Given limitations of available analytical tools for exact derivation of limit cycles in discrete-time models, we approximate the time spent in the ghost attractor τ by considering a power law for the time spent in a limit cycle (Medeiros et al. 2017): Whenever n M > n ∧ and p M < p ∧ , System 5 shows that n k+1 > n k and p k+1 < p k for all k > M. Therefore, we identify the time the system escapes the ghost attractor In this figure we consider p 0 = 10, n 0 = 1, δ p = 0.9, γ p = 0.9912γ * p , σ = 2.67, δ n = 0.8, γ n = 1, R = 2. Although we know that the system will converge to the equilibrium point, this convergence takes over 5000 time steps. as τ = min{M : n M > n ∧ , p M < p ∧ }. Figure 6 shows that this approximation using the power law provides a reasonable approximation.

Discussion
In this work we have identified two types of long transients, crawl-by transients and ghost attractors, that can appear in a consumer-resource system with group defense with discrete reproductive pulses. Our long-term dynamics are qualitatively different from those found in Cui et al. (2016), where they identified a variety of bifurcations and chaotic dynamics. The key differences between the two models are that, while that of Cui et al. (2016) models reproduction as a continuous process and integrates the density-dependent growth for a case with overcompensation in discrete time, our model considers reproduction as a discrete event and has a saturating density-dependent function. The model in Cui et al. (2016) is a discrete-time model similar to the Ricker When modeling reproduction as a discrete process with saturating (Beverton-Holt style) density dependence rather than overcompensation, our analysis did not suggest that increased reproduction numbers leads to instabilities in our model.
In addition, discrete reproduction events are one of the main reasons we see the long transients analyzed in this model. The crawl-by transient observed at high consumer densities (Theorem 1) is caused by a sudden crash in the adults of the resource population, which is followed by a slow crash of the consumer population due to its inability to find enough resource for self-replacement. Although the adult resource population is almost nonexistent, the few remaining individuals eventually lead to an increase the resource population when the consumer population becomes small enough.
The other reason long transients appear in this model is due to the self-replacement of consumers depending on the ability of resource to defend themselves. The Type IV Holling functional response produces a bifurcation on the proportionality constant γ p at the value γ * p given by Eq. 6. This constant can be associated with the conversion capability of consumers, i.e., the amount of energy invested in reproduction activities. When the conversion capability of consumers is too small (γ p < γ * p ), group defense of resource will prevent self-replacement of consumers at high resource densities, which will lead to collapse of consumers. When this conversion capability is high enough, self-replacement can be satisfied, and the consumer-resource cycles of Fig. 3 will occur. These cycles and their condition for existence resemble those found in other models where a mechanism of group defense of resource is considered (Ajraldi et al. 2011;Venturino 2011;Venturino and Petrovskii 2013).
In the case where the system presents consumer-resource cycles, the resourcedominated phase will include a crawl-by transient when the conversion capability is close to this bifurcation value (Theorem 3). This will follow by a crash of the resource population, where the consumer-dominated phase appears and presents the crawl-by transient previously described. This behavior presents an alternative perspective to the concept of alternate stable states (Beisner et al. 2003), where the different "alternative stable states" constitute long transients, which may resemble stable states during a long period of time, which then transition into the other phase and stay in a different long transient. In reality, stochasticity may render this juvenile survival when rare impossible or accelerate consumer death, which may lead the model to a stable state in a shorter period of time (Reimer et al. 2021).
When the conversion capability of consumers is smaller than the critical value γ * p but close to it, these quasiperiodic orbits do not disappear completely and stay as ghost attractors. This ghost attractor stays until the resource density surpasses a given threshold (the equilibrium value n ∧ ) and the system enters the basin of attraction of the resource-only equilibrium. The emergence of this ghost attractor is caused by group defense, because in its absence (σ = 0), the unstable equilibria that cause the quasiperiodic orbits ( p ∨∧ , n ∨∧ ) do not exist. In their absence, resource population density will consistently increase and the consumer density decrease. The estimation of the transient time of the ghost attractor shows one of the limitations of our analysis, as the theory to study limit cycles in discrete-time systems is not developed enough to precisely analyze the transient time of this ghost attractor. Given the seasonality of the reproduction and recruitment for many organisms (Arreguin-Sanchez 1992;Cameron 1986;Russell et al. 1977;Wallace 1985), a continuous-time model may not properly reflect the biological dynamics we are interested in. Despite this challenge, the expression for the transient time found for transient limit cycles in continuous-time systems in Medeiros et al. (2017) is a reasonably accurate fit in our model. The transient times of the ghost attractor found in our work are similar to those found in the predator-prey model with group defense of Venturino and Petrovskii (2013). However, our biological mechanisms differ, as their transients could be attributed to search time of prey from the predators through space, a feature not explicitly modeled in our work. In contrast, the length of the ghost attractor in our model can be attributed to the length of the crawl-by transients that are part of the cycle itself, which are periods of low population growth for either the consumer or the resource.
In conclusion, we show how long transients can appear in predator-prey systems with group defense and discrete recruitment pulses. A possible extension of this model is to explicitly consider the dynamics of the juvenile stages through a continuous-time model, which could give a more accurate approximation of the transient times found in this paper. A multi-stage model would also allow exploration of the effect of relative adult versus juvenile vulnerability to consumption on the transient dynamics.
by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

Appendix: A Fixed points of System 5 and their stability
The fixed points of System 5 ( p, n) satisfy the equations p = δ p p + γ p pn 1 + σ n 2 n = δ n n exp − γ n p 1 + σ n 2 + Rn exp(− p) 1 + n .
If p = 0, then the second equation gives us two solutions for n, n = 0 and If p = 0, then the first equation has two solutions for n given by where n ∨ corresponds to the solution with a − sign and n ∧ to the solution with a + sign. These solutions are positive whenever γ p ≥ 2 √ σ (1−δ p ). In such case, plugging n ∨∧ into the second equation provides us with the following expression: Then there is an unique value p ∨∧ that solves the trascendental equation This equation in p has an unique solution as the function is monotonic for p and satisfies lim p→−∞ f ( p) < 0 and lim p→∞ f ( p) > 0. For it to be biologically relevant, we also require lim p→0 f ( p) < 0, which will occur when or, after reorganizing the terms, n ∨∧ < n * . The Jacobian of the system J is the following: From here, the extinction equilibrium satisfies which has eigenvalues δ p < 1 and δ n + R > 1. Therefore the extinction equilibrium is a saddle. For the resource-only equilibrium, the upper right term of the Jacobian matrix equals 0 whenever p = 0. Therefore J (0, n * ) is a triangular matrix, with the eigenvalues being the diagonal terms Because R > 1 − δ n , 0 < λ 2 < 1. λ 1 , on the other hand, λ 1 will produce a change in stability when In this case, the equilibrium is stable whenever γ p < γ * p and a saddle when γ p > γ * p . Because n * > 1/ √ 2σ , plugging γ p = γ * p in Eq. 12, we have that n * = n ∧ . Based on the conditions for ( p ∧ , n ∧ ) to be biologically reasonable, this implies that at γ p = γ * p , the system goes through a transcritical bifurcation, where the carrying capacity (0, n * ) changes stability.
The trascendental equation that describes p ∨∧ renders it impossible to analyze them directly. However, a numerical exploration in Fig. A.6 shows that these equilibria are unstable for γ p < γ * p and ( p ∧ , n ∧ ) becomes stable for γ p > γ * p . Combining this result with the condition for the equilibrium point ( p ∧ , n ∧ ) to be biologically relevant (Eq. 15), we find that when γ p > γ * p , there are no stable fixed points in the first quadrant (i.e., R 2 + ). This transcritical bifurcation occurs with almost any combination of parameters in our region of interest. To show this, we perform a similar analysis as those in Khan et al. (2016); Murakami (2007). When γ p = γ * p , we can rewrite our system in diagonal form and centered around the origin by making the change of variables: Fig. 7 Numerical values of different equilibria for n and their stability as we vary γ p in our region of interest. The red points correspond to unstable equilibria, whereas the blue points correspond to stable equilibria. In this figure, δ p = 0.9, σ = 2.67, δ n = 0.8, γ n = 1, R = 2.
Plugging o(ε) into the formula for p 2 gives us that: In addition, if n m = o(ε), the equation for n m+1 satisfies: n m+1 n m = O (δ n exp (−γ n p m ) + R exp(− p m )) .
While p m = O(ε −1 ), this expression will satisfy n m+1 /n m < 1. We can thus assume that the expression is satisfied until p m+1 = O(1). Therefore, when p m = O(ε −1 ), the consumer population time evolution can be approximately solved as This expression stops working when p m = O(1), and thus resource will start a recovery afterwards. We can estimate the order of magnitude of such m by plugging p m = 1 above. Solving for m gives us that m = log(ε) log(δ p ) which is the expression that proves the theorem.

Appendix: C Proof of Theorem 2
To show the existence of this theorem, we use Theorem 2.9 of Magal and Zhao (2005).
To do this, we consider the first quadrant M as a metric subspace of R 2 with metric d(x, y) = x − y 2 induced by the Euclidean norm. Let T : M → M be given by We show that T is a point dissipative, compact map on M. Because M is a subspace of R 2 , compactedness is trivial. A map is point dissipative if there is a bounded set and eigenvectors This system has for solution the expression where a, b are constants. If we let m = 0, then we can solve the linear system which has solutions In particular, this gives us that a = ε u . Because γ p > γ * p , then λ 1 > 1, and λ 2 < 1. Therefore, for big m, System 40 can be approximated as Thus, System 5 will stay near the resource-only equilibrium as long as (x m , y m ) = O(ε). In particular, the System will escape the saddle point when x m = O(1). Plugging in x m = 1 into the approximated solution lets us find M that solves the equation This has for solution M = log 1 ε log(λ 1 ) which is the expression that proves the theorem.