Locust Dynamics: Behavioral Phase Change and Swarming

Locusts exhibit two interconvertible behavioral phases, solitarious and gregarious. While solitarious individuals are repelled from other locusts, gregarious insects are attracted to conspecifics and can form large aggregations such as marching hopper bands. Numerous biological experiments at the individual level have shown how crowding biases conversion towards the gregarious form. To understand the formation of marching locust hopper bands, we study phase change at the collective level, and in a quantitative framework. Specifically, we construct a partial integrodifferential equation model incorporating the interplay between phase change and spatial movement at the individual level in order to predict the dynamics of hopper band formation at the population level. Stability analysis of our model reveals conditions for an outbreak, characterized by a large scale transition to the gregarious phase. A model reduction enables quantification of the temporal dynamics of each phase, of the proportion of the population that will eventually gregarize, and of the time scale for this to occur. Numerical simulations provide descriptions of the aggregation's structure and reveal transiently traveling clumps of gregarious insects. Our predictions of aggregation and mass gregarization suggest several possible future biological experiments.

These equations generalize the classic swarming model which describes a single population density field advected by a velocity field arising from social interactions.Eq. ( 2) has been studied extensively in one and two spatial dimensions for various social interaction functions represented by Q, whose negative gradient is the effective social force [?, ?, ?, ?]. Depending on Q, solutions include steady swarms, spreading populations, and contracting groups (i.e., blow-up) [?,?,?].
In our two-phase model Eqs.(1), the velocities are and the social interaction potentials Q s,g are Here, R s , R g , A g are interaction magnitudes and r s , r g and a g are interaction length scales.We require R g a g −A g r g > 0 and A g a 2 g −R g r 2 g > 0 so that Q g includes short range repulsion and long range attraction, as in [?, ?, ?], as this is the clumping regime, appropriate to capture the tendency of gregarious locusts to aggregate.We model the density-dependent rates of interconversion of the solitary and gregarious forms as The parameters δ 1,2 are maximal rates and k 1,2 are characteristic locust densities at which the transitions occur at half of their maximal values.To the best of our knowledge, our work is the first to consider locust phase changes via continuum modeling of locust density [?, ?, ?, ?].

Parameter selection and estimation
As discussed in the main text, for our numerical results, we use two different sets of phase change parameters.For both sets, we use the same social interactions parameters, and we now describe our choices for these.To estimate R s , R g , and A g , we use explicit velocity computations.The speed of a locust when it is alone varies between 72 -216 m/hr, while the speed of a locust in a group varies in a tighter range of 144 -216 m/hr [?].To make a rough estimate of R s , we imagine a hypothetical semi-infinite density field ρ(x, y) = ρ group H(x) where H(x) is the Heaviside function and, as mentioned in the main text, ρ group = 65 locusts/m 2 is the approximate critical density of a gregarious group [?].A solitary locust placed at the swarm's edge (at the origin) should move to the left with maximal velocity v max s = −216 m/hr.From Eqn. (3), which we solve to find R s = 11.87 m 3 /(hr • locust).Similarly, a gregarious locust at the origin should move to the right with maximal velocity v max g = 216 m/hr, so A gregarious locust placed to the left of the swarm at a distance equal to the attraction length scale a g = 0.14 m should also move to the right, but with a slower velocity which we take to be the minimal velocity in a crowd, v min g = 144 m/hr.Thus v g (−0.14, 0) = {−∇Q g * ρ group H(x)} (−0.14,0) = v min g .
These two conditions determine R g = 5.13 m 3 /(hr • locust) and A g = 13.33 m 3 /(hr • locust) In the main text, we present numerical simulations of Eqs.(1) in one spatial dimension.For these simulations, we take δ 1,2 , r s , r g , and a g as above, since these parameters do not depend on spatial dimension.For the remaining parameters, we follow a process similar to that described above, and choose k 1,2 = k = 8 locusts/m, R s = 6.83 m 2 /(hr • locust), R g = 6.04 m 2 /(hr • locust), and A g = 12.9 m 2 /(hr • locust).

Homogeneous steady states
For any set of initial conditions, the mean locust density ρ 0 is known, and corresponds to the total density at the homogeneous steady state (HSS).Accordingly, there is a family of homogeneous steady states parameterized by ρ 0 .The corresponding solitary and gregarious HSS components, obtained by setting time and space derivatives to zero in Eqs.(1) are When we later consider stability of homogeneous steady states, it will be convenient to discuss the fractions φ s,g of solitarious and gregarious locusts, where φ s + φ g = 1.Using Eqn.(9), we know that for homogeneous steady states, Here, γ = δ 1 /δ 2 is the ratio of maximal solitarization rate to maximal gregarization rate, K = k 1 /k 2 is the ratio of the characteristic solitarization and gregarization densities for individuals, and ψ = ρ 0 /k 2 is a rescaled density.Note that φ g is monotonically increasing in ψ, and hence in ρ 0 ; that is to say, as total density increases, the gregarious fraction increases.

Linear stability analysis
To study the stability of the HSS in Eqs.(9), we consider small perturbations s 1 , g 1 about s 0 , g 0 s(x, t) = s 0 + s 1 (x, t), g(x, t) = g 0 + g 1 (x, t), (11) so that ρ(x, t) = s 0 + g 0 + s 1 (x, t) + g 1 (x, t).Substituting Eqn.(11) into Eqn.(1) and expanding to first order in the perturbations, we find the linearized equations where Here, A, B > 0 for all ρ 0 > 0 since f 1 is a monotonically increasing function of ρ 0 and f 2 is a monotonically decreasing one.To further analyze the linearized equations, we Fourier expand the perturbations as We allow for an infinitely large domain so that there are no restrictions on q; in other situations, q must be suitably restricted in order to satisfy boundary conditions.Substituting Eqn.(14) into Eqn.( 12) yields ordinary differential equations for each Fourier mode amplitude.We write these in matrix form, Here, q = |q| is the perturbation wavenumber, and Q s,g (q) are the Fourier transforms of the two dimensional social interaction potentials, The eigenvalues λ 1,2 (q) of L(q) are Since λ 2 < 0, instability occurs only when λ 1 > 0. For convenience, we rewrite λ 1 in terms of the gregarious mass fraction φ g , Now we factor out the attractive part of the gregarious term, namely This yields Since the prefactor is negative, and we seek conditions for a positive eigenvalue (signifying growth of perturbations, and hence instability), we focus on when the term in square brackets becomes negative.The dependence on φ g occurs via the prefactor (1 − φ g )/φ g in front of a positive term.For possible instability, this term should be small, meaning that φ g should be sufficiently large (since this prefactor is monotonically decreasing with φ g ).Since φ g increases monotonically with ρ 0 (as discussed above), instability may occur as ρ 0 is increased.We now show that instability first occurs at the wavenumber q = 0 (meaning that perturbations that first lead to instability are long wavelength).We again focus on the bracketed quantity in Eq. (21).If this term becomes negative, it must do so for the value of q at which the first two terms are (together) minimized, since these are positive terms and the negative term, −1, is a constant.It is biologically reasonable to assume that a g ≥ r s (with equality achieved for our chosen social interaction parameters).Therefore, the first term is either constant or monotonically increasing in q.It is also biologically reasonable to assume that a g > r g , in which case the second term is monotonically increasing in q.Thus, the first two terms together are monotonically increasing in q, so their minimum occurs at q = 0, and this will be the first wavenumber to trigger instability.Thus, if we are looking for the instability that occurs as φ g increases, it is sufficient to consider what happens at q = 0.
We substitute q = 0 into the bracketed term in Eqn.(21) and ask for what value of φ g the resultant expression changes sign (to find the threshold level of gregarious locust fraction needed for instability).Setting that bracketed term to zero we obtain Instability is achieved for values of φ g greater than this threshold value.
To obtain a more explicit condition for instability in terms of the density ρ 0 , we substitute φ * g into Eq.(10), which relates gregarious fraction to total (scaled) density.Rearranging, we obtain the biquadratic equation where For any biologically meaningful solutions, the solution for ψ 2 must be positive.From the quadratic formula, we have Since A > 0 and C < 0, the discriminant is positive.Hence, for the plus sign choice, ψ 2 > 0. For the minus sign choice, ψ 2 < 0 and hence we eliminate this possibility.The final result for the critical scaled density is This is the result that we use to produce instability contours in the K-γ plane (Fig. 2 in the main paper).

Numerical simulation method
We simulate Eqs. ( 1)-( 5) in one spatial dimension.We use periodic boundary conditions on a domain of length L with a fine grid consisting of N = 1024 points (necessary to resolve the steep edges of clusters that form).To approximate an unbounded domain, one may take the limit of large L. The social interactions Q s,g in ( 4) must be adapted to be commensurate with a periodic domain.We begin with the function Q(x) = e −|x|/r , which is the building block of Q s,g .We calculate the discrete Fourier transform F of −∂ x Q on our domain as where r is the decay length scale in Q and ∆ = L/N is the grid spacing.From Eqn. (27) it is straightforward to compute the Fourier transforms of Q s,g .Convolutions are equivalent to products in Fourier space, providing excellent computational savings (and thus justifying the choice of a periodic domain).We compute velocities by convoluting the density with −∂ x Q s,g pseudospectrally.The flux term in Eqs. ( 1) is instead evaluated via a fourth-order accurate central finite difference.
The emergence of discontinuities in s and g causes ringing in the pseudospectral evaluation of the velocity term.In order to smooth this effect, we incorporate small amounts of numerical diffusion.Another standard approach would be to incorporate high wave number filtering in the simulation.We choose numerical diffusion because it also serves as the macroscopic description of random motion, which locusts certainly display.We implement diffusion in a split-step manner, alternating with the dynamics of Eqs.(1)- (5).Time-stepping is performed with the fourth-order Runge-Kutta method.We also threshold our velocity field at every time step so that it does not exceed v max g .Without this thresholding, individual locusts achieve velocities of up to approximately 1.5 times v max g at an intermediate stage of our simulation.It is crucial to point out that this thresholding only affects the speed of the transient clumps; it does not affect the initial instability (which is small amplitude, and thus has a small velocity) and similarly, it does not affect the late-stage bulk dynamics (which are nearly spatially stationary).

Introduction
Outbreaks of locusts such as Schistocerca gregaria, Locusta migratoria, and Chortoceites terminifera regularly afflict vast areas of Northern Africa, the Middle East, Asia, and Australia.Depending on climate and vegetation conditions, billions of voracious locusts aggregate into destructive swarms that span areas up to a thousand square kilometers.A flying locust swarm can travel a few hundred kilometers per day, stripping most of the vegetation in its path [1][2][3][4].A recent locust plague in West Africa (2003)(2004)(2005) severely disrupted agriculture, destroying $2.5 billion in crops destined for both subsistence and export.Despite control efforts totalling $400 million, loss rates exceeded 50% in certain regions [5,6].These numbers alone attest to the urgency of finding better ways to predict, manage, and control locust outbreaks.
Between outbreaks, locusts are mainly antisocial creatures who live in arid regions, laying eggs in breeding grounds lush with vegetation.Resource abundance may, on occasion, support numerous hatchings, leading to a high population density.Overcrowding at resource sites promotes transition to a social state in a self-reinforcing process.The social locust nymphs may display mass migration behavior.Within the newly formed group, individuals cohere via sensory communication, whether visual, chemical, and/or mechanical [3].Outbreaks may be exacerbated in periods of drought, when large numbers of locusts congregate on the same breeding or feeding grounds [7][8][9].
Locusts are phase polyphenic: while sharing the same genotype, individuals may display different phenotypes [10,11] that incorporate variations in morphology [12], coloration [13], reproductive features [14] and, significantly, behavior [15,16].An individual can change from a solitarious state (preferring isolation) to a gregarious one (seeking conspecifics).Behavioral state is plastic [3,11,15] and strongly dependent on local population density: in sparse surroundings, a gregarious locust transitions to the solitarious state [15] and vice versa in crowded environments.These phase transitions are called solitarization and gregarization.Gregarization dominates when large numbers of locusts gather at the same site, potentially leading to a destructive outbreak [8,9].
Locust gregarization may be induced by visual, olfactory, or tactile cues.For the desert locust Schistocerca gregaria, the most potent stimulus is tactile: repetitive stroking of the femora of hind legs [15][16][17] functions as a crowding indicator.Mechanosensory stimulation of leg nerves leads to serotonin cascades in the metathoracic ganglion, and initiates gregarious behavior [16][17][18].Gregarization can be induced by rubbing a locust's hind leg for 5 s per minute during a period of 4 hr [17].Cessation of physical contact leads to solitarization after 4 hr, though the degree of solitarization achieved during that time depends on the individual's ancestry.
Experiments and models have shed much light on how group alignment [19][20][21][22] and group motion [23,24] depend on group size or density and treatments such as diet and denervation.For instance, a lowprotein diet (which motivates cannibalism in locusts) leads to stronger interactions between individuals and lowers the threshold density beyond which mean speed and group coherence increase [24].Other data-driven studies include models based on a well-known physics paradigm for self-propelled particles [25] and explore the transition between a disordered and a coherent marching group.Both [26] and [27] study the dynamics of rolling patterns formed by flying, gregarious swarms.A logistic map was introduced in [28] to describe phase change via a birth rate and a carrying capacity dependent on population density modulated by stochastic effects.
Our current work complements previous locust modeling studies in several ways.First, many of the previous models are individual-based (Lagrangian) simulations, where the position, velocity, and interactions of individual locusts are tracked [19,20,[22][23][24].Ours is density-based (Eulerian), allowing techniques of partial differential equations (PDEs) and their extensions (integro-PDEs) to be utilized.Second, we concentrate on gregarious-solitarious transitions not yet explicitly considered in [20,24].We address intrinsic attractive-repulsive social interactions, whereas many current models consider interactions with clumped resources and environmental heterogeneity as their focal points [8,9].Finally, some models [24] include anisotropic interactions such as different responses to anterior and posterior neighbors, or consider Newtonian dynamics.To explore minimal mechanisms sufficient for band formation, our work instead uses isotropic interactions and a kinematic approach.The open problem we address via mathematical modeling is to quantify and describe collective gregarization, a key, early process that necessarily occurs before the emergence of a destructive locust outbreak.We do this by linking the physiology of individual-level phase change and interphase interactions to predictions at the level of the gregarious hopper band as a whole.
We investigate the onset of an outbreak by constructing a continuum mathematical model of behavioral phase for interacting gregarious and solitarious locusts.We classify and quantify group dynamics in wide swaths of parameter space, a task which is challenging by numerical techniques alone.We find that in the limit of low densities, both phases are uniformly spread and the solitarious phase dominates.For sufficiently large populations, a dense, traveling patch of gregarious locusts suddenly emerges, while solitarious locusts become more and more scarce.We identify locust clustering at high densities with the onset of a hopper band.Through analysis of our model, we calculate the critical density beyond which the gregarious group forms, and for the final ratio of gregarious to solitarious locusts.We determine these quantities in terms of behavioral parameters at the level of individual locusts, hence connecting individual and group properties.Our model also displays population-level hysteresis, which has implications for locust management.

Model Model construction
Locusts in a group are subject to attractive and/or repulsive forces based on combined sensory, chemical, and mechanical cues that affect their motion.We assume that sensing is directionally isotropic, a reasonable approximation [29] for organisms receiving sensory inputs of a variety of types, although directional models are possible as well [30].Rather than tracking individual locusts, we consider a population density field ρ(x, t) moving at velocity v(x, t).Continuum population modeling [31,32] allows us to apply analytical tools in order to characterize swarm formation and structure.Our work draws from classic swarm modeling in which a conserved population density field ρ moves at a velocity v that arises from social interactions: This is the well-known mass balance equation that tracks individuals moving collectively at velocity v.It is typically assumed that individuals can sense the population density nearby, and that this sensing gives rise to attractive-repulsive social forces F, or alternatively, social potentials Q (the negative gradients of which are forces).Within this context, the contribution ρ(x ′ , t) of a small clump of individuals at location x ′ to the force on the individual at position x is given by F The corresponding velocity is proportional to the forces exerted by neighbors at all spatial locations, so that v(x, t) is given by integration over all The expression for the velocity v(x, t) in Eqn. ( 2) is a convolution of the density ρ(x, t) and the social interaction force −∇Q(x − x ′ ), which describes the influence of the locust population at location x ′ on that at location x.This is a common formulation of so-called nonlocal interaction models [33][34][35][36], which capture interactions that are spatially distributed, in contrast to pure partial differential equations, which include only local terms such as derivatives and gradients, and which describe interactions only over infinitesimal ranges.Nonlocal aggregation models have been studied for various social interactions Q; known solutions include steady swarms, spreading populations, and contracting groups.We use the notation v = −∇Q * ρ to denote the convolution in Eqn.(2).We assume that Q(x − x ′ ) is radially symmetric and depends only on the distance between x and x ′ .The detailed forms of Q in the case of solitarious and gregarious locusts will be described later.
To adapt Eqs. ( 1) and ( 2) to biphasic insects, we introduce separate density fields for solitarious and gregarious locusts, s(x, t) and g(x, t), respectively, and the total local density ρ = s + g.With marching locusts in mind, we consider a two-dimensional geometry, with Ω representing the spatial domain and x = (x, y) as spatial coordinates.We now include the phase transitions between solitarious and gregarious locusts.To do so, we define two density-dependent functions, f 1 (ρ) for the the rate of gregarious-tosolitarious transition, and f 2 (ρ) for the rate of solitarious-to-gregarious transition.Our model thus reads where the velocities are given by These equations are complete once we specify the solitarious and gregarious social interactions Q s,g and the density-dependent conversion rates f 1,2 .Since solitarious locusts are crowd-avoiding, we take Q s to be purely repulsive.Gregarious locusts, on the other hand, are attracted to others, except for short-distance repulsion due to excluded volume effects.Hence, we model Q s and Q g as where R s , R g , A g are interaction amplitudes that determine the strengths of attraction and repulsion, and r s , r g and a g are interaction length scales that represent typical distances over which one locust can sense and respond to another.The above forms of Q s,g describe social interactions that decay exponentially away with distance from the sensing individual and are chosen to be isotropic for simplicity.As evident from Eqn. (5), Q s is purely repulsive for all choices of R s and r s .On the other hand, Q g is the difference of two exponentials, implying that there may be a distance at which repulsion and attraction balance, resulting in no net contribution to the velocity.The location of this balance point can be obtained by imposing −∇Q g = 0 to obtain the critical distance Depending on the choice of social interaction parameters, the expression for d may yield unphysical results such as negative distances.The distance d also pertains only to two isolated locations x and x ′ and does not capture population-level features.Even for meaningful values of d, a collection of individuals interacting under Q g may disperse, aggregate, or clump.It is thus important to choose the appropriate parameter ranges for a g , r g , A g and R g so that the tendency of gregarious locusts to aggregate is modeled properly.Mathematical studies have shown that in order for cohesiveness to occur, the parameters in Q g must lie in a particular regime that leads to clumping [37].Thus, we require R g a g − A g r g > 0 so that repulsion dominates at short length scales, and A g a 2 g − R g r 2 g > 0 so that attraction dominates at longer ones.Taken together, these conditions guarantee a meaningful critical distance d and macroscopic clumping behavior.We assume these conditions to hold for the remainder of this paper.
It remains to specify how density affects transitions from one phase to another.We call upon the biological observation that at higher densities, gregarization proceeds more quickly and solitarization more slowly.We model the phase conversion rates with the rational functions The parameters δ 1,2 are maximal phase transition rates and k 1,2 are characteristic locust densities at which f 1,2 take on half of their maximal values.Note that f 1 decreases with ρ, capturing the inverse relationship between solitarization rate and density, while f 2 increases with ρ and saturates at δ 2 , describing speedier gregarization at higher densities.
Our complete model consists of Eqs. ( 3)-( 7) together with initial conditions specifying s(x, 0) and g(x, 0).We consider a spatially periodic domain, which simplifies both numerical simulation and mathematical analysis.In certain laboratory studies using ring-shaped arenas, such boundaries are natural (while being less ideal for comparison with field studies) [20].We do not include locust reproduction or death as these occur on much longer time scales than phase change.
The model presented here is a general one containing some fundamental elements of locust dynamics.This work can be readily modified and extended to include details pertaining to different locust species, interactions with the surrounding environment, locust reproduction, and more.For instance, in our model, we have not explicitly accounted for the differing activity levels of solitarious and gregarious individuals [11].Additionally, while gregarization is relatively fast for Schistocerca gregaria, full solitarization may occur only after several generations of locusts.The phase conversions of Chortoicetes terminifera, on the other hand, are characterized by similar timescales for the two phase conversions, so that both gregarization and solitarization occur rapidly within the lifetime of a single locust individual [38].On another note, vegetation or waterway patterns may impose spatial inhomogeneities such as non-uniform initial distributions of solitarious locusts, or attraction to preferred sites.Preexisting models in the literature have pointed out the important link between the spatial distribution of vegetation, as well as nutritional quality, on locust clustering, gregarization, and swarming [8,9,39,40].All of these elements could be used to refine our model for predictive purposes.However, as the first work in the continuum modeling of locust population phase change, ours begins with the fundamental model contained in Eqs. ( 3)- (7).Our model is complementary to the preexisting ones in that we focus on how inherent inter-individual interactions can lead to gregarization and swarming, even in a spatially homogeneous environment.Multi-generational dynamics, differential activity levels, resource distribution, and related factors could be considered as possible extensions of our work.

Parameter selection
Some of our results are analytical formulas, which may be evaluated for any desired parameter values.Other results depend on numerical computations, and these require specific choices of parameters.For these results, we consider two different sets of phase transition parameters.(1) Most of our numerical results have been obtained using our default set of parameters, based on estimates from the biological literature.Specifically we take δ 1,2 = δ = 0.25 hr −1 , corresponding to a gregarization time scale of approximately 1/δ = 4 hr for desert locusts (for whom some -but not total -solitarization occurs on the same time scale) [11,17].We also take k 1,2 = k = 65 locusts/m 2 , since for desert locusts, the critical density for the onset of collective motion is 50 -80 locusts/m 2 [24].We will allow for some deviation from δ 1 = δ 2 and k 1 = k 2 via a parameter sensitivity analysis.(2) To examine situations with large differences in the rates of gregarization and solitarization, we consider an alternative set of parameters with δ 1 = 0.025 hr −1 and δ 2 = 0.25 hr −1 , so that gregarization is an order of magnitude faster that solitarization.We take k 1 = 20 locusts/m 2 and k 2 = 65 locusts/m 2 to model a gregarious-to-solitarious transition that occurs at a higher density threshold than the solitarious-to-gregarious transition.
We use the same social interaction parameters for all results (variations from this set are accounted for by a sensitivity analysis).To estimate the social interaction length scale parameters in Eqs.(5), we apply the results of [20,24], which identify the "sensing range" of a desert locust as 0.14 m, and the "repulsion range" as 0.04 m, close to the approximately 0.05 m body length of a desert locust at the fifth instar of its development.For the gregarious phase we thus set the repulsion length scale at r g = 0.04 m and the attractive one at a g = 0.14 m, corresponding to the experimental sensing range.These choices agree with theoretical studies showing that for cohesive swarms, attraction occurs over longer length scales than repulsion [33,41].We also assume that solitarious locusts are repelled from others at their sensing range, so that r s = 0.14 m.These choices satisfy r g < a g = r s which is assumed for the remainder of this paper.
Finally, we estimate R s , R g , and A g via explicit velocity computations.The speed of a locust when it is alone varies between 72 -216 m/hr, depending on diet [24].At the upper end, this is roughly one body length per second.When it is moving in a group, the individual's speed varies in a tighter range of 144 -216 m/hr [24].In making our phase-dependent velocity estimates, we interpreted the "moving alone" and "moving in a group" data as typical to solitarious and gregarious locusts, respectively.Using these biological measurements and Eqn.(4), we find R s = 11.87 m 3 /(hr • locust), R g = 5.13 m 3 /(hr • locust), and A g = 13.33 m 3 /(hr•locust).Details are given in Text S1.Our choices of social interaction parameters satisfy conditions mentioned in the previous section, namely R g a g − A g r g > 0, and A g a 2 g − R g r 2 g > 0 so that gregarious insects will clump.
Most of our parameter choices have been inferred or estimated from published laboratory experiments.It is possible however, that in the field, some parameter values may be quite different from the ones we have used.For instance, locusts in the field may pause while marching to perch on the vegetation, giving rise to an effective speed that is lower than what measured in lab experiments, where perching does not occur.It is also noteworthy that gregarious locusts are more active than solitarious locusts, a fact that is reflected by our method of choosing R s , R g , A g from estimates of the velocities of individuals when moving alone and in a group.As we describe below, we analyze our model varying all parameters within reasonable bounds: our results are qualitatively the same.

Results
We first determine the simplest solutions to the model, namely those for which the densities of gregarious and solitarious locusts are in a spatially uniform steady state.We probe the stability of that uniform state using linear stability analysis (LSA), a calculation that addresses whether small, spatially nonuniform perturbations grow or decay.This is equivalent to determining the signs of eigenvalues of the linearized system, where positive (negative) eigenvalues imply growing (decaying) perturbations.The rate of initial growth/decay depends on the wavenumber of the perturbation.The growing perturbations can be interpreted in terms of nascent aggregates of locusts, and the wave numbers as the number of aggregates per unit area.The analysis provides a condition for the onset of aggregation, namely the emergence of positive eigenvalues of the linearized model.In our case, this aggregation condition is shown below in Eqn.(14).LSA cannot, in general, predict the ensuing dynamics once perturbations have grown to a large size.Further analysis uses an approximation to eliminate the spatial dependence of the model, which enables an analytical prediction of the proportion of solitarious and gregarious locusts on a longer time scale.To visualize the dynamics of aggregation, we perform numerical simulations in one spatial dimension using the linear stability analysis to identify regimes of interesting behavior.The model displays population-level hysteresis.

Homogeneous steady states
The solitarious s 0 and gregarious g 0 homogeneous steady-state (HSS) solutions of Eqn.(3) can be written in terms of the total uniform density ρ 0 , which is simply the mean value of ρ for a specified initial condition.The full expressions for s 0 and g 0 in terms of ρ 0 appear in Text S1; in the small ρ 0 limit these are approximately while in the limit of large ρ 0 we find The low density HSS is thus composed mostly of solitarious locusts and vice versa for the high-density case, showing the non-monotonicity of s 0 with respect to total density ρ 0 .In Fig. 1(A) we plot the HSS s 0 (middle solid blue curve) and g 0 (middle broken green curve) for our default set of phase change parameters, k = 65 locusts/m 2 and δ = 0.25 hr −1 .As shown, s 0 initially increases with ρ 0 .At a critical density ρ * , s 0 reaches a maximum, whereas g 0 keeps increasing monotonically.Fig. 1(B) shows a blow-up of the region near ρ * .For our default parameters, the maximum value s max 0 is attained at ρ * = k, the same density value for which solitarious and gregarious densities coincide so that s max 0 = s 0 (ρ * ) = g 0 (ρ * ) = k/2.However, this feature is a result of our choice k 1 = k 2 and δ 1 = δ 2 .In general, the point of maximum solitarious density and the point of equal solitarious and gregarious density do not coincide, as is directly deducible from the full expressions for s 0 and g 0 in Text S1.To give a sense of detuning from our parameter estimates, we also calculate and plot s 0 and g 0 for parameter sets chosen randomly from uniform distributions centered at our estimated default set of values for {δ 1 , δ 2 , k 1 , k 2 }.The bottom and top curve in each set show the 25th and 75th percentile values.
We also study a much more general case where δ 1 = δ 2 , k 1 = k 2 , in keeping with the distinct rates of transition and critical transition densities seen biologically.As an alternative way to understand the HSS solutions, we consider the fractions φ s,g of solitarious and gregarious locusts, where φ s + φ g = 1.As shown in Text S1, for the HSS, Here, γ = δ 1 /δ 2 is the ratio of maximal solitarization rate to maximal gregarization rate, K = k 1 /k 2 is the ratio of the characteristic solitarization and gregarization densities for individuals, and ψ = ρ 0 /k 2 is a rescaled spatially homogeneous density.The gregarious fraction φ g is monotonically increasing in ψ, and hence in ρ 0 ; that is to say, as total density increases, the gregarious fraction increases.For small ρ 0 , φ g ≈ 0, but as ρ 0 increases, there is a crossover between solitarious and gregarious populations.Uniformly spread solitarious populations cannot be sustained when the density is too high: the gregarious state will necessarily become the dominant one.

Linear stability analysis
To determine conditions under which a nearly uniformly spread locust population aggregates or disperses, we study the linear stability of the HSS (details appear in Text S1).The calculation is a standard but somewhat tedious exercise.In nonlocal systems such as ours, linear stability results depend on the Fourier transforms Q s,g (q) of the interaction potentials Q s,g .For our locust model, the stability of the HSS depends on the eigenvalue where q = |q| is the perturbation wave number and the Fourier transforms Q s,g (q) in two dimensions are Observe that the eigenvalue λ 1 (q) depends on all of the individual-based parameters governing rates of phase change (via s 0 and g 0 ) and all of the social interaction amplitudes and length sensing length scales.
The HSS derived in the previous section is stable to small perturbations if λ 1 (q) < 0 for all q.If λ 1 (q) > 0 for some q, then the HSS is unstable to perturbations of those wave numbers.
Our full analysis of this eigenvalue appears in Text S1.We formulate the instability condition in terms of φ g , If this condition is satisfied, initially small perturbations from the uniform steady state will grow.This inequality is a key result, and implies that if a sufficiently large fraction of the population is gregarious, the HSS solution is unstable.To obtain a more explicit condition in terms of the density ρ 0 , one must substitute φ * g into Eqn.(10), which relates gregarious fraction to total (scaled) density.One may then calculate the critical density ρ 0 above which the HSS is unstable.Since φ g and ρ 0 are monotonically related, we conclude that the HSS solution is unstable for sufficiently dense populations.The algebra is tedious, and relegated to Text S1.Instead, we present a contour plot in Fig. 2 which succinctly illustrates the stability features of the HSS.The phase change parameter ratios γ = δ 1 /δ 2 and K = k 1 /k 2 vary along the horizontal and vertical axes and the contours indicate the critical value of rescaled density ψ * = ρ * 0 /k 2 .For scaled densities greater than ψ * , the HSS solution is unstable.The critical scaled density is monotonically increasing in both γ and K. (Note that for an accurate biological interpretation, one must multiply ψ * by k 2 in order to obtain the unscaled critical density ρ * 0 .)Upon inserting our default parameters in Eqn. ( 14) we find that the homogeneous solution is unstable for ρ 0 > ρ * 0 = 62.3 locusts/m 2 .This value corresponds to the left border of the grey region in Fig. 1.For ρ 0 > ρ * 0 , to the right of the border, we expect the onset of a locust hopper band, i.e., formation of patches of high locust density that can seed the clustering and gregarization of other locusts.In Fig. 1, linear instability can occur even at densities ρ 0 for which s 0 exceeds g 0 for our chosen parameters (represented by the center solid blue and center broken green curves).This result implies that the onset of instability leading to mass gregarization can take place even if solitarious locusts initially outnumber gregarious ones.We will later discuss mass gregarization in more detail.To visualize detuning from this set of parameters, we include the 25th and 75th percentile values of ρ * 0 for onset of instability as vertical purple lines; these are again calculated by drawing 10,000 random samples of the parameters k 1,2 , δ 1,2 , R s,g , r s,g , A g , and a g .As seen from Fig. 1(b) our conclusions are robust across the randomly chosen parameter sets.
For our default set of biological parameters, φ * g ≈ 0.479 via Eqn.( 14) and ρ * 0 turns out to be near k = k 1,2 .We stress that generically, it is not the case that ρ * 0 needs to be near k 1 and/or k 2 .For our default parameter set, K = γ = 1, in which case ψ = 0.959, so that the critical value ρ * 0 is 95.9% of k 2 , namely 62.3 locusts/m 2 .However, for different choices of K and γ, drastically different outcomes are possible.For instance, for our alternative parameter set where K ≈ 1/3 and γ = 1/10, the critical density is ρ * 0 = 15.9 locusts/m 2 , which is quite disparate from the individual gregarization density of 65 locusts/m 2 , and is also less than the solitarization density of 20 locusts/m 2 .Furthermore, for different choices of the social interaction parameters entering into Eqn.(14), it is possible to obtain a critical gregarious fraction φ * g that is much less than 1/2, meaning that instability and clumping can occur even with just a few gregarious insects.
For ρ 0 > ρ * 0 , we can also find the wave number q max corresponding to the most rapidly growing perturbation.Fig. 3 shows q max for our chosen parameters (center curve) as well as the 25th and 75th percentile values over the 10,000 random parameter draws.The most unstable wave number q max grows rapidly as a function of ρ 0 and then saturates at q max ≈ 8.89 m −1 , corresponding to a length scale 2π/q max ≈ 0.71 m and indicating that the most quickly growing perturbations occur on the length scale of a few locust bodies.
Our linear stability analysis describes the behavior of small perturbations of uniform steady states, and is not expected to predict long-term or large-amplitude dynamics.For large perturbations, linear analysis is void.Additionally, even to analyze small perturbations of states other than uniform steady states, a different analysis would be needed.

Numerical simulation
To illustrate the swarm dynamics described by Eqn.(3), we simulate the model on a one-dimensional periodic domain of length L = 3 m for a total population of M = 50 locusts.Periodicity of the domain is an important aspect of a robust numerical platform devised for these simulations: we exploit the fact that convolutions Q * ρ are easy to compute in Fourier space (where they are simply products, i.e., Q• ρ), which significantly reduces the computational overhead.Computational issues associated with such convolutions also restrict us to one-dimensional simulations at present.At t = 0 all locusts are solitarious and are randomly perturbed from the uniform density s = M/L, where M is the total population mass We adjust some parameters so as to adapt our model to the one-dimensional case.Specifically, one must take square roots of k 1,2 in order to collapse densities in a square to densities along a line segment.Consequently, for our default parameter set we choose k 1,2 = k = 8 locusts/m and δ 1,2 = δ = 0.25 hr −1 , whereas for the alternative set we use k 1 = 4.5 locusts/m, k 2 = 8 locusts/m, δ 1 = 0.025 hr −1 and δ 2 = 0.25 hr −1 .In both cases we take the interaction amplitudes R s = 6.83 m 2 /(hr • locust), R g = 6.04 m 2 /(hr • locust), and A g = 12.9 m 2 /(hr • locust), which have also been adapted from their original values to the one-dimensional case.The interaction length scales r s , r g , and a g are the same as for the two-dimensional case.Details of the numerical method and the parameter choices appear in Text S1.
Results are shown in Fig. 4 for the default parameter set and in Fig. 5 for the alternative set.In each case, the snapshots show s(x, t) (dashed blue curve) and g(x, t) (solid green curve) at selected times.Starting from the randomized solitarious state at t = 0 hr, locusts rapidly redistribute to a roughly spatially uniform density until t ≈ 3 hr.Tiny variations are present but not visible on the scales of these figures.Gregarization and subsequent rapid spatial segregation follow.In Fig. 4, between t ≈ 3.42 hr and t ≈ 3.47 hr, two compactly supported clumps of gregarious locusts emerge, superposed on a background of sparse, solitarious individuals.A similar transition occurs between t ≈ 3.15 hr and t ≈ 3.17 hr in Fig. 5, but for these parameter values, we find initial clustering with three, rather than two density peaks.The number (or alternatively, length scale) of transient clumps that form appears to be selected dynamically.This intermediate dynamical selection process and the coarsening that ensues are avenues for future numerical and analytical investigation.In each example, the disjoint clusters quickly merge due to the long-range attraction of gregarious individuals.A single remaining pulse is formed by t ≈ 3.49 hr in both cases and travels until t ≈ 6.5 hr, at which time the majority, but not all, of the solitarious locusts have transitioned to the gregarious form.Gregarization continues during the subsequent hours, albeit at a slower rate.For both figures, the gregarization of the final clump continues slowly, approaching an equilibrium at exponentially long times.
To study the locust gregarization process further, we define the total mass of solitarious and gregarious locusts, S and G, as so that the total population mass is M = S + G.We also define the mass fractions which we before calculated for HSS solutions, but we now generalize for spatially varying states.These quantities will be useful to further our mathematical analysis.Fig. 6 shows φ s (t) (blue curve) and φ g (t) (green curve) as arising from the numerical simulations depicted in Fig. 4 and Fig. 5. Several distinct regimes are visible, and we discuss these below.

Spatially-homogeneous and spatially-segregated bulk theories
As visible in the second and third panels of Fig. 4 and Fig. 5, the early-time dynamics of Eqs. ( 3) are approximately spatially homogeneous.As a result, spatially-dependent terms in Eqs.(3) are negligible, ρ is approximately constant, and hence the governing equations are linear ordinary differential equations (ODEs) that are easily solved.We write the solution of these ODEs in terms of the mass fractions φ s,g , where we have used the initial condition φ s (t = 0) = 1.This analytical solution is plotted in Fig. 6 as a dotted line, and agrees closely with the numerical results for the first few hours.
In the later panels of Fig. 4 and Fig. 5, gregarious and solitarious locusts spatially segregate into areas with disjoint support.This means that in each distinct region, ρ(x, t) ≈ s(x, t) or ρ(x, t) ≈ g(x, t).We thus consider a bulk model reduction to study the dynamics of the two non-overlapping solitarious and gregarious populations.In particular, we assume that solitarious locusts are spread throughout most of the domain Ω, covering an area denoted α s , whereas gregarious locusts are confined to a region with area α g .Within these areas, local densities are approximately S/α s and g = G/α g .By integrating Eqs.(3) over the domain and assuming that s and g are approximately constant in their support, we obtain where The numerical solution of these ODEs (dashed lines in Fig. 6) agrees closely with the late time fullscale numerical simulation results, where we use values of α s,g measured empirically from the terminal equilibrium.One can reduce Eqn.(19) to a single nonlinear ODE using φ s = 1 − φ g , though this equation is not amenable to analytical solution.Since we are interested in the large population limit for which we expect potential large scale gregarization, we instead study Eqs.(19) for large M .In this case, to leading order in M , the bulk model reduces to Given the expressions for c 3,4 and the fact that φ s , φ g ≤ 1, the first term is O(1) whereas the second one is much smaller, O(1/M 2 ).For large M then, and to leading order, φ s decays exponentially in time with rate δ 2 .This result is based on the assumption of a segregated state, and thus would be expected to occur only once segregation is nearly complete.
Since for large M (nearly) the entire population will eventually become gregarious, the critical density ρ * 0 is a crucial result.If the population is in the stable regime (where ρ 0 < ρ * 0 ) then mass gregarization can be avoided and solitarious and gregarious locusts can coexist as uniformly spread populations.However, as soon as the population shifts beyond the border of stability (where ρ 0 > ρ * 0 ) the group gregarizes and the onset of a locust hopper band is inevitable.

Phase change and hysteresis
The biological literature discusses the importance of hysteresis in locust phase change, as reviewed, for instance, in [11].It is important to disambiguate the possible meanings and interpretations of phase change hysteresis, to place this phenomenon within the context of our model, and most especially, to distinguish between hysteretic features at the individual and population levels.
One type of hysteresis is simply defined as "rates of gregarization [that] differ from rates of solitarization" [11].Within our model, this type of hysteresis may be interpreted as cases where δ 1 = δ 2 or k 1 = k 2 .Our results thus far have accounted for this type of hysteresis in three ways.First, for our primary parameter set in which δ 1 = δ 2 and k 1 = k 2 , we have allowed deviations from equality by performing a sensitivity analysis incorporating variations of up to 30% from the base parameter values, as represented in the results of Fig. 1 and Fig. 3. Second, for our alternative parameter set, we have chosen δ 2 = 10δ 1 and k 2 ≈ 3k 1 .And finally, for analytical results such as the homogeneous steady states and their stability, we have obtained analytical formulas into which any values of δ 1,2 and k 1,2 can be substituted.
Another interpretation of hysteresis relates to "solitarization [having] two phases: an initial rapid phase and a second, slower phase that requires insects to be maintained in isolation across successive moults -or generations" [11].Our model is constructed on the time-scale of a single generation, and thus we cannot account for this type of hysteresis, which would require a multi-generational model.
Finally, we can consider population-level hysteresis.In the context of our model, this type of hysteresis refers to macroscopic properties of solutions of Eqn.(3) (which are outputs of the model) as opposed to differences in individual-level parameters (which are inputs to the model) as in the first type of hysteresis described above.Numerical results suggest that our model has population-level hysteresis; see Fig. 7.This figure shows the gregarious mass fraction φ g as the average density ρ 0 (total mass M divided by domain length L) is varied as a control parameter.All phase change, social interaction, and physical domain parameters are the same as in Fig. 5.
The solid (dashed) red curve is an analytical result, representing the stable (unstable) HSS solution, as calculated previously via linear stability analysis.For small values of ρ 0 , the HSS is stable to small perturbations.If locusts join the initially stable population the average density ρ 0 will increase (assuming a fixed spatial domain), shifting the uniform state to the right along the red curve; as yet no clustering will be evident.Beyond the point labeled with an asterisk, the uniform HSS loses stability and clustering occurs, as previously described.This corresponds to a jump represented by the vertical black arrow.The clustered state (green) is now stable.We next ask what happens if locusts are now removed from the aggregate, which corresponds to a reduction in ρ 0 (moving to the left in Fig. 7).We answer this question numerically, by gradually subtracting mass from the population, allowing the system dynamics to evolve, and plotting the gregarious fraction as a function of mass.As the mass is slowly removed, the solution tracks leftwards along the green curve, indicating the persistence of the gregarious band.In fact, the band persists even partway into the regime where the HSS is linearly stable.
This dynamically observed hysteresis suggest that (for our model) a gregarious aggregation cannot be eliminated by reducing overall density to a low enough level where the HSS is linearly stable.This result has implications for locust control, as we discuss below.

Discussion
In this paper, we derived, analyzed, and simulated a model for the movement, social interactions, and density-dependent interconversions of the solitarious and gregarious forms of phase polyphenic locusts.The model is based on experimental observations and measurements, parameter values inferred from preexisting work, and basic assumptions about individuals' rules of behavior.We included social exchanges via repulsive and/or attractive interactions for gregarious and solitarious individuals, and we accounted for phase change with density-dependent transitions, with crowding favoring solitarious-to-gregarious conversions.Our model was formulated in terms of continuum equations, allowing us to apply classical techniques such as linear stability analysis and bulk approximation.Since these methods were applied in two spatial dimensions, our results are relevant to insects aggregating in two dimensional structures such as hopper bands.We also provided example simulations in one spatial dimension as proof of principle, and as an indication of typical dynamics.
Our model explicitly takes into account intrinsic social interactions between individuals, in contrast to pre-existing models that focus on how insects respond to quality and spatial heterogeneity of nutrition or other environmental factors [8,9,24].These approaches are complementary, showing that both intrinsic and extrinsic factors that affect local densities also affect the gregarization transition.
Many of our results are achieved via mathematical analysis.The power of mathematical analysis is that it creates an explicit connection between individual-level and group-level quantities, e.g., via the inequality Eqn.(14).Once we identify the sensing range and interaction strength parameters in Eqn.(5) which govern individual locust attraction to and repulsion from others, we are able to calculate the critical density beyond which mass gregarization occurs.
Briefly, our results and predictions can be summarized as follows: (1) Locusts exist in a spatially uniform steady state distribution only up to a critical total population density.(2) Beyond this critical density, the uniform distribution can not be maintained, and massively dense gregarious clusters form.
(3) Linear stability analysis allows us to understand how the critical density depends on dimensionless ratios of the biological parameters.This dependence is summarized in Fig. 2. Our analysis also yields the most unstable cluster spacing (from the wave number of the most unstable modes).( 4) Numerical simulations illustrate the rapid transitions that take place once gregarization is initiated.Dense packs of gregarious locusts form and grow, and these move and sweep up solitarious locusts in their vicinity.
(5) Via bulk approximation, we find estimates for the long-time mass fraction dynamics of solitarious and gregarious locusts.In the large population limit, the entire population will become gregarious.Bulk theory and simulations agree well, as shown in Fig. 6. (6) Our model displays population-level hysteresis, via which the critical density at which a gregarious aggregation forms from a dispersed population can be significantly higher than the density at which a gregarious aggregation would break up, as shown in Fig. 7.
Our results shed light on locust control strategies in two ways.First, given the mass gregarization that takes place past the point of linear instability, the density threshold for this instability is a crucial quantity.In accordance with the idea proposed in [42], our work identifies a threshold below which populations should be kept in order to avoid a gregarious outbreak (assuming biological parameters are known to a sufficiently accurate degree).Furthermore, we have shown how this population-level property depends on individual-level parameters, finding a nontrivial relationship.Second, the apparent population-level hysteresis shows that dispersing a gregarized band, perhaps by killing individuals with pesticides, is harder than preventing group formation in the first place in that band annihilation requires a significantly lower locust density.In short, hysteresis implies that prevention could be more easily achieved than control.
Like all models, ours has its limitations.We did not include features of the environment such as vegetation, shown to have important influence on local crowding and hence gregarization.Our simplifications lead to mathematical tractability, while limiting the direct biological relevance of the model at present.In the field, locusts encounter patchy vegetation and other environmental influences, and adding such factors to the model would make it more relevant to field experiments.Since we have not explicitly included resource gradients or other environmental cues, we do not here recapitulate the long-range motion of locust bands, but merely their formation and clustering.Including environmental factors constitutes an extension of the current framework.Similarly, simulations in two spatial dimensions are more challenging and remain open for future investigation.
Our work suggests several future biological experiments.First, as always, more accurate knowledge of model inputs would lead to better results.For our model, key inputs include the social interaction parameters, namely the length scales (r s , r g , and a g ) and interaction amplitudes (R s , R g , and A g ) in Eqn.(5) that we inferred from careful experiments such as those in [24].However, to our knowledge, most of these parameters have not been directly measured in experiments on individuals.Second, we encourage observations of macroscopic group properties that could be compared to outputs of our model.These outputs include densities and sizes of bands.Additional quantitative field measurements along the lines of [43] could help validate and refine our model.Finally, we can imagine experiments that would probe important aspects of the system dynamics (as opposed to physical properties of the bands themselves).Hopper bands are known to undergo complicated dynamics, including splitting and merging [3].BBC video shows an example of such phenomena in Locustana pardalina bands [44].More accurate data for the dynamics of wild groups, including times for group formation and distances between merging bands and tributaries, could be compared to clumping time and length scales identified by our model.We are especially curious about experiments in which the critical average density for population-level gregarization and clumping might be probed in a controlled lab experiment, perhaps by slowly adding solitarious individuals into a large arena.Experimental measurements like those we have mentioned here would also motivate future two-dimensional extensions of our model where the streaming dynamics of hopper bands, the effects of the environment, and other stimuli could be more fully explored.(A) Spatially homogeneous solitarious (s 0 , solid blue) and gregarious (g 0 , broken green) steady state locust density as the total locust density (ρ 0 ) is varied.For each set of curves, the middle curve represents the solution for our default phase change parameters, k 1 = k 2 = 65 locusts/m 2 and δ 1 = δ 2 = 0.25 hr −1 .The bottom and top curve in each set show parameter sensitivity; they are the 25th and 75th percentile values for s 0 and g 0 over 10,000 parameter sets of {δ 1 , δ 2 , k 1 , k 2 } sampled from uniform distributions centered at the default values and varying by ±30%.In both the thin grey and red regions, the HSS is linearly unstable to small perturbations.Additionally, in the red region, g 0 > s 0 , while in the grey and white regions, the opposite holds.(B) A blow-up of the boxed transition region in (A) around which the value of g 0 overtakes s 0 .The dashed black vertical lines indicate the 25th and 75th percentile for this transition.The solid purple vertical lines indicate the 25th and 75th percentile values for the onset of linear instability.For rescaled densities greater than that value, the HSS solution is unstable.The critical rescaled density is monotonically increasing in both γ and K.The arrow along the horizontal axis indicates the direction γ moves if the relative rate of gregarization is increased (faster gregarization).The arrow along the vertical axis indicates the direction K moves if the relative density threshold for gregarization is decreased (easier gregarization).For an accurate biological interpretation, one must multiply ψ * by k 2 in order to obtain the unscaled critical density ρ * 0 . .Maximally unstable perturbation wave number q max for homogeneous steady states with total density ρ 0 .to Fig. 1, the middle, bottom, and top curves show results for the 25th and 75th percentile as computed from 10,000 random parameter draws centered around our default parameter set.At low densities, there are no unstable perturbation wave numbers.Just past the critical density ρ * 0 , q max increases rapidly and then plateaus.For our default parameters, q max asymptotes to 8.89 m −1 corresponding to a length scale of 0.71 m.  (3), as can be seen up until the second row of panels, where the small instability becomes significant.Two compactly supported clumps of gregarious locusts form, superposed on a very sparse population of solitarious insects.In the third row, the gregarious group travels as a propagating pulse, and eventually stops.During this stage, the gregarious and solitarious populations are essentially non-overlapping in space.As shown in Fig. 6, the group continues to slowly gregarize after it becomes stationary.Population-level hysteresis as a function of average density ρ 0 .Gregarious mass fraction φ g as the average density ρ 0 (total mass M divided by domain length L) is varied as a control parameter.We use our alternative set of phase change and social interaction parameters, as in Fig. 5.The solid (dotted) red curve represents the stable (unstable) homogeneous steady state solution, as calculated via linear stability analysis.As ρ 0 passes through the point of linear instability (marked with an asterisk) the solution jumps up to the green curve, which represents compactly supported gregarious aggregations obtained via numerical simulation, similar to the final states of Fig. 4 and Fig. 5.As ρ 0 is decreased by slowly subtracting mass from aggregations on the green curve, the system remains on the upper branch even for values of ρ 0 sufficiently small as to be in the regime where the uniform state is stable, thus demonstrating dynamical population-level hysteresis.

Figure 1 .
Figure1.Spatially homogeneous steady states (HSS).(A) Spatially homogeneous solitarious (s 0 , solid blue) and gregarious (g 0 , broken green) steady state locust density as the total locust density (ρ 0 ) is varied.For each set of curves, the middle curve represents the solution for our default phase change parameters, k 1 = k 2 = 65 locusts/m 2 and δ 1 = δ 2 = 0.25 hr −1 .The bottom and top curve in each set show parameter sensitivity; they are the 25th and 75th percentile values for s 0 and g 0 over 10,000 parameter sets of {δ 1 , δ 2 , k 1 , k 2 } sampled from uniform distributions centered at the default values and varying by ±30%.In both the thin grey and red regions, the HSS is linearly unstable to small perturbations.Additionally, in the red region, g 0 > s 0 , while in the grey and white regions, the opposite holds.(B) A blow-up of the boxed transition region in (A) around which the value of g 0 overtakes s 0 .The dashed black vertical lines indicate the 25th and 75th percentile for this transition.The solid purple vertical lines indicate the 25th and 75th percentile values for the onset of linear instability.

Figure 2 .
Figure 2. Linear stability of spatially homogeneous steady state (HSS) solutions.The dimensionless phase change parameter ratios γ = δ 1 /δ 2 and K = k 1 /k 2 vary along the horizontal and vertical (we have used log axes).The contours indicate the critical value of rescaled density ψ* = ρ * 0 /k 2 .For rescaled densities greater than that value, the HSS solution is unstable.The critical rescaled density is monotonically increasing in both γ and K.The arrow along the horizontal axis indicates the direction γ moves if the relative rate of gregarization is increased (faster gregarization).The arrow along the vertical axis indicates the direction K moves if the relative density threshold for gregarization is decreased (easier gregarization).For an accurate biological interpretation, one must multiply ψ * by k 2 in order to obtain the unscaled critical density ρ * 0 .

1 )Figure 3
Figure 3. Maximally unstable perturbation wave number q max for homogeneous steady states with total density ρ 0 .to Fig.1, the middle, bottom, and top curves show results for the 25th and 75th percentile as computed from 10,000 random parameter draws centered around our default parameter set.At low densities, there are no unstable perturbation wave numbers.Just past the critical density ρ * 0 , q max increases rapidly and then plateaus.For our default parameters, q max asymptotes to 8.89 m −1 corresponding to a length scale of 0.71 m.

Figure 4 .
Figure 4. Numerical simulations of Eqs.(3).Snapshots depict the numerical solution of Eqs.(3)-(7) at different times t (in hr) on a periodic domain of length L = 3 m with the default set of phase change parameters.See also Fig.5for a comparison with the alternative parameter set.The solitarious (gregarious) density (in locusts/m) as a function of spatial position (in m) is shown in blue (green).The total population mass is M = 50 locusts and the initial condition is set at g(x, t = 0) = 0 and s(x, t = 0) given by a random perturbation centered around s = M/L.The top row of panels shows the fast smoothing of the initial state, and the subsequent evolution.Gregarization (approximately) occurs according to the spatially homogeneous version of Eqn.(3), as can be seen up until the second row of panels, where the small instability becomes significant.Two compactly supported clumps of gregarious locusts form, superposed on a very sparse population of solitarious insects.In the third row, the gregarious group travels as a propagating pulse, and eventually stops.During this stage, the gregarious and solitarious populations are essentially non-overlapping in space.As shown in Fig.6, the group continues to slowly gregarize after it becomes stationary.

Figure 5 .Figure 6 .
Figure 5. Numerical simulations of Eqs.(3).Similar to Fig.4, snapshots at different times t (in hours), but for the alternative set of phase change parameters.Note that three, rather than two clumps of gregarious locusts form at intermediate times.This simulation is continued until t = 80 hr (last frame) to show the stability of the final cluster of gregarious locusts.

Figure 7 .
Figure 7. Population-level hysteresis as a function of average density ρ 0 .Gregarious mass fraction φ g as the average density ρ 0 (total mass M divided by domain length L) is varied as a control parameter.We use our alternative set of phase change and social interaction parameters, as in Fig.5.The solid (dotted) red curve represents the stable (unstable) homogeneous steady state solution, as calculated via linear stability analysis.As ρ 0 passes through the point of linear instability (marked with an asterisk) the solution jumps up to the green curve, which represents compactly supported gregarious aggregations obtained via numerical simulation, similar to the final states of Fig.4and Fig.5.As ρ 0 is decreased by slowly subtracting mass from aggregations on the green curve, the system remains on the upper branch even for values of ρ 0 sufficiently small as to be in the regime where the uniform state is stable, thus demonstrating dynamical population-level hysteresis.