Agent-based and continuous models of hopper bands for the Australian plague locust: How resource consumption mediates pulse formation and geometry

Locusts are significant agricultural pests. Under favorable environmental conditions flightless juveniles may overpopulate and aggregate into coherent, aligned swarms referred to as hopper bands. These bands are often observed as a propagating wave having a dense front with rapidly decreasing density in the wake. A tantalizing and common observation is that these fronts slow and steepen in the presence of lush, edible vegetation. This suggests the collective motion of the band is mediated by resource consumption. Our goal is to model and quantify this effect. We focus on the Australian plague locust, for which excellent field and experimental data is available. Exploiting the alignment of locusts in hopper bands, we concentrate solely on the density variation perpendicular to the front. We develop two models in tandem; an agent-based model that tracks the position of individuals and a partial differential equation model that describes locust density. In both these models, locust are either stationary (and feeding) or moving. Resources decrease with foraging. The rate at which locusts transition between moving and stationary (and vice versa) is enhanced (diminished) by resources abundance. This effect proves essential to the formation, shape, and speed of locust hopper bands in our models. From the biological literature we estimate ranges for the ten input parameters of our models. Sobol sensitivity analysis yields insight into how the band's macroscopic characteristics vary with changes in the input parameters. By examining 4.4 million parameter combinations, we identify biologically consistent parameters that reproduce field observations. We thus demonstrate that resource-dependent behavior can explain the density distribution observed in locust hopper bands. This work suggests that foraging should be an intrinsic part of future modeling efforts.


Introduction
Locusts are a significant agricultural pest in parts of Africa, Asia, Central and South America, and Australia. They form large groups of millions to billions of individuals that move collectively, rapidly consuming large quantities of vegetation [1,2]. Their collective movement occurs during both nymphal and adult stages of development and is associated with an epigenetic phase change from a solitary to a gregarious social state which is mediated by conspecific density and abiotic factors [1,[3][4][5][6]. In the earlier stages of collective motion, flightless nymphs, also called hoppers, form aligned bands that march along the ground, often through agricultural systems where they cause significant crop damage as they forage and advance [4,7,8]. Some species, such as the brown locust Locustana pardalina, form intertwining streams of relatively homogeneous density [1,2,8]. By contrast the Australian plague locusts Chortoicetes terminifera form wide, crescent-shaped bands that contain a high density in front and a rapidly decreasing density behind [4,9,10]. As bands of C. terminifera move through a field, they create a sharp transition from undamaged vegetation in front of the band to significant defoliation immediately behind the band (see Fig 1). In natural systems, they tend to consume one of several species of grasses; in agricultural systems, they tend to eat primarily pasture and sometimes early stage winter cereals [11]. Our goal in this study is to model the formation of hopper bands for the Australian plague locusts and to argue that the interaction between feeding and motion is likely responsible for their characteristic form.  [4,9,10], with a steep leading edge and shallower decay behind that is roughly exponentially decreasing in density [9]. Notable in aerial photographs [10] of these hopper bands propagating through crops is the contrast between the verdant green of the unperturbed crops in front of the band and the lifeless brown in the pulse's wake. The one meter wide strip above represents the dimensions we use to model locust movement in a single dimension, as described below.
Endemic to Australia, the Australian plague locust C. terminifera is the most common locust species on the continent. For ease, we will henceforth refer to C. terminifera simply as "locust". Locust development relies on rainfall events which result in the rapid growth of green plant biomass needed for locust survival [12]. If enough successive rainfall events occur followed by a lack of rain, locusts migrate and breed in such a way that high egg density occurs in a single geographic location leading to outbreaks of emergent nymphs [13]. At sufficient densities, these locusts undergo an epigenetic phase change from a solitary state (in which locusts actively avoid conspecifics) to a gregarious state characterized by attraction to conspecifics, leading to aggregation and collective motion [3]. Nymphs then form aligned marching groups referred to as hopper bands, which may move tens to hundreds of meters per day [4]. While marching, they consume large quantities of green biomass with an individual eating approximately one third to one half of its body weight per day [11]. Locusts proceed through five nymphal stages, called instars, with marching behavior beginning during the second instar [4]. Approximately four weeks after eggs hatch, the locusts reach adulthood and are then capable of forming highly destructive mobile flying swarms [14]. To study marching, we focus on the third and fourth instar.
Marching behavior is strongly mediated by ambient temperature [10]. When temperatures are too low, such as during the early morning, individuals stop to bask in the sun. When temperatures are in an optimal range, bands become marching occurs. During the hottest part of the day, soil temperatures can become too high and the locusts may stop to roost on cooler vegetation [10]. During periods of sustained marching, when locusts at the front of a hopper band encounter food resources, they stop and often feed [4,10]. Locusts farther back in the band may continue to move forward, eventually passing those that stopped. This creates a "leap-frogging" type motion with a cycling of individuals in the dense front of the band. Clark [4] describes this exact behavior in his 1949 paper: Those hoppers behind the front were in places which had been partly or wholly eaten out, and thus lacked the same stimulus of food to stop them. As their average rate of progress was greater than that of the hoppers in the front, they tended to overtake them, becoming in turn slowed down in their progress by the presence of food.
Thus, locust motion in the band relies on resource density and individuals stopping to forage; however, we are aware of no models of Australian plague locust band movement that have incorporated foraging behavior. This motion is consistent with the concept of the Marginal Value Theorem [15] which quantifies optimal foraging: if a cost is assigned to foraging that is proportional to the resource density, it costs less energy per unit resource to move when sufficiently higher density resources are available.
Previous modeling efforts have considered both agent-based and continuous models [1]. The majority of these models have focused on social behavior, notably alignment, attraction and repulsion with respect to conspecifics [16][17][18][19][20][21][22]. Many of the agent-based models considered the pause/go behavior of locusts [19,21,22], and other insects [23]. Continuous models have been used to model transition between stationary and moving states [24] and gregarization [25]. Foraging has been modeled as an agent-based model [26] and resource distribution effects on peak density has been modeled as an energy minimization problem [27].
By contrast, our models concentrate on the interaction of locust pause and go motion with food resources. We assume that hoppers march in an aligned band through a field of finite resources, which is depleted as the locusts stop to feed. We then develop a model for the probability of movement as a function of resource availability, considering in tandem an agent-based model (ABM) which tracks individual locusts and a partial differential equation (PDE) model that considers locust densities. Both models produce traveling pulse solutions that are consistent with the detailed field observations of Buhl et al [9]. The ABM is easily simulated, allows us to track individuals within the swarm, and captures the natural stochasticity of a biological process. By contrast, the PDE produces smooth solutions and lends itself to analysis and a detailed characterization of how observable outcomes such as pulse geometry, pulse speed, and density of resources left unconsumed in the pulses wake are related to the model's parameters.
Our main conclusion will be that foraging and resource mediation of stationary-moving transition rates produces the pulse structures observed in the field supporting the observations and inferences of Clark [4]. We note that this is in contrast to previous studies [17,22] who suggest that the formation of sharp asymmetric fronts can be explained solely by social forces. A strength of our model is that it reproduces the observed profiles [9] quantitatively from biologically realistic parameters.
In the Section 1 we construct our two models beginning with biological and simplifying assumptions, and ending with parameter identification from empirical field data in Table 2. Section 2 contains our results, namely that both models produce a traveling pulse of locusts precisely when the locusts' stationary-moving switching is dependent upon the amount of nearby resources. Evidence consists of numerical simulations for the ABM, traveling wave analysis for the PDE, and a robust sensitivity analysis of the models to changes in the input parameters. In Section 3 we revisit our main findings and outline extensions of this work incorporating more biological complexity. Our appendices in the Supplemental Material contain mathematical analysis and proofs substantiating results for the PDE.
1 Models and Methods

Basic Assumptions
We demonstrate that the formation of a locust hopper band traveling through a landscape may be driven by resource consumption. We construct two minimal models to investigate this relationship and, in doing so, make a series of assumptions simplifying the biological reality.
• We assume hopper bands consist of flightless nymphs and behaviorally identical in all regards.
• During a typical day, a hopper band has periods of collective movement (marching) and collective resting (roosting when temperatures are too high or too low) [1,4,8,9]. We model only the periods of marching.
• We assume that while the group is marching, it is made up of individuals in either a stationary state or a moving state. Further, only stationary locusts feed while moving locusts propagate forward with a constant speed that represents an average of crawling and hopping.
• We assume that locusts transition between stationary and moving states at a rate depending solely on the resources nearby. These transition processes are completely memoryless, which implies that locusts experience neither hunger nor satiation.
• The transition rate from stationary to moving is positive and decreases as the resource density increases. Similarly, the transition rate from moving to stationary is positive and increases with increasing resource density. Below we model these rates as decreasing/increasing exponential functions, which allows their values to saturate at large resource concentrations (see Fig 2).
• Locusts are known to align their direction of movement with their nearest neighbors and may align with environmental cues such as wind or the location of the sun [2,4,8]. Therefore, we assume individuals move parallel to one another, creating a constant direction of movement for the entire band. This alignment is the only social force we incorporate. In all other respects, we assume locusts act independently.
• Since there is no lateral movement, we model behavior in a narrow strip aligned to the direction of movement, as shown in Fig 1. For dimensional consistency of the model, we assume the transverse width to be 1 meter wide.
• Since locusts feed much more quickly than vegetation grows, we assume that resources (food) can only decrease. Moreover, resources are identical so that they can be characterized by a single variable. Prior to locust arrival, we assume available resources are spatially constant.
It is notable that in this model there is no explicit social interaction between the locusts; interaction is mediated solely through the consumption of resources. Social interaction is included implicitly in that our one-dimensional model assumes that all the locusts are aligned with a common direction of motion. Our primary goal with this study is to show that linking resource consumption and pause-and-go motion can explain the observed hopper band morphologies. Adding social interactions to this model might improve its efficacy in one dimension and would be an absolute necessity for modeling behavior in two-dimensions.

General Model Formulation
Within the framework described above, we build two models: an agent-based model (ABM) which tracks individual locusts and a partial differential equation model (PDE) that determines locust density. These models share much in their basic structure. Table  1 compares their independent and state variables and Table 2 lists their common model parameters.
In the ABM, space and time lie on a discrete, evenly spaced lattice (x n , t m ) while in the PDE space and time (x, t) are continuous. In both models, S and M denote the number or density of stationary and moving locusts respectively. For the ABM, the number of stationary (moving) locusts at x n ,t m is denoted S n,m (M n,m ). For the PDE, the analogous continuous quantities for the density of locusts are S(x, t) and M (x, t).
Resources edible by locusts are measured by the non-negative scalar density variable R; specifically, the resource density in the agent-based model is R n,m , and the resource density in the continuous model is R(x, t).
We assume that the group rate of feeding is proportional to the product of the stationary locust density and the resource density; that is, (rate of change of resources at a given location) = −λSR (1) where λ is a positive rate constant that describes how quickly individual locusts consume resources. This implies that an individual locust's foraging efficiency decreases as resources become scarcer at their location. This is not an explicit implementation of the Marginal Value Theorem but fits the general concept of foraging efficiency within a patch decreasing due to searching time, not satiation by the forager [15]: as the resources at a location are eaten, locusts have difficulty locating the next unit to consume, reducing the overall rate of resource consumption at that location. For the rates of transition between moving and stationary, we use common exponentially saturating functions of resources as illustrated in Fig 2. The stationary to moving rate is denoted k sm while the moving to stationary rate is called k ms . Specifically, where γ, δ > 0, 0 ≤ β ≤ θ, and 0 < η ≤ α. The conditions on the parameters guarantee that k sm (R) is a decreasing function and k ms (R) is a increasing function of R. (Most of the analytical results concerning the PDE model hold for any choice of monotone switching rates -see Supplemental S2 Appendix for details.) This functional form derives from the assumption that the transition rate's sensitivity to changes in resources is proportional to the resource availability [15]. Transition rates for stationary to moving k sm (gold) and moving to stationary k ms (purple) with α = 0.12, β = 0.02, γ = 0.03, δ = 0.015, η = 0.005, and θ = 0.14.

Agent-Based Units Continuous Units Description
Model Model Biologically, this implies that when encountering excess resources, there will be a high proportion of stationary locusts, and doubling the excess resources will do little to change the proportion of stationary locusts. Similarly, when resources are scarce, locusts are most likely to transition from stationary to moving and least likely to stop. Mathematically, this functional form preserves the positivity of the transition rates and means that the transition rates are constant in the limit of abundant resources. In the PDE model, the transition rates k sm , k ms appear as coefficients in growth and decay terms in the differential equations. In the ABM we use a stochastic version of these transitions. At each time step, locusts switch from stationary to moving via a transition probability p sm and from moving to stationary via p ms , both of which are functions of R n,m . The smooth transition rates k ms and k sm can be understood to be derived from these probabilities as the time step ∆t approaches zero. Working backward and assuming ∆t is small, we will utilize the following approximations, This is equivalent to assuming that each locust undergoes only a single transition in any given time step.

Agent-Based Model (ABM): Pause and Go Motion on a Space-time Grid
We now construct our agent-based model (ABM) describing the behavior of each individual locust. The temporal evolution of the ABM may be thought of as a probabilistic cellular automaton. The model is one dimensional in space, representing a 1-meter-wide cross section of the locust hopper band. Our ABM tracks the position of each locust, their states (stationary or moving), and the spatial availability of resources (food). Locust position and the spatial distribution of resources are confined to a discrete lattice of points given by x n = n∆x and time t m = m∆t, for n, m ∈ N. We fix ∆x = v∆t so that a moving locust moves forward one step on the lattice each time step.
Let X i (t m ) be the position of the i th locust at time t m . Let σ i (t m ) be a binary state variable where σ i = 1 when the locust is moving and 0 otherwise. The motion of the locusts can now be expressed succinctly as where we have applied the value of the state variable at t m throughout the interval of length ∆t. Note this artifice ensures that the values X i remain on the lattice for all time t m . We model transitions between stationary and moving states with a discrete-time Markov process given via the probabilities in Eq (3). Thus, at time t m , each locust at position x n has a p sm (R n,m ) probability to switch from stationary σ i = 0 to moving σ i = 1 or a p ms (R n,m ) probability to switch from moving to stationary.
To consider the macroscopic behavior, we define histogram variables by counting the number of locusts with each state at each space-time grid point. Let S n,m = X n (t m )(1 − σ n (t m )) = # of stationary (σ i = 0) locusts at (x n , t m ). (6) We model the resources with a scalar variable R n,m which is defined as available food, measured in grams, at time t m in the interval of width ∆x centered at x n . Following Eq (1) and converting S n,m to a density, we have Solving Eq (7) (assuming S n,m is constant between t m and t m+1 ) yields R n,m+1 = R n,m e −λSn,m ∆t ∆x .
Biologically, this evolution implies that the resources in a patch of vegetation infested by a group of stationary, feeding locusts will decrease by approximately half in an amount of time inversely proportional to λ times the number of locusts in the patch. That is, the half-life of resources in the patch is ln(2)∆x/(λ · # locusts). We initialize each simulation with R n,1 = R + , indicating an initially constant field of resources. Together with initial conditions, Eqns. (3), (4), and (8) specify the evolution completely. Our agent-based model then takes the form of three sequential, repeating steps for each locust agent: 1. Update state S or M according to the Markov process.
2. If in state M , move to the right ∆x.
3. If in state S, decrease resources in current location.
Each locust performs each of these steps simultaneously with all other locusts, and resources in each location are also updated simultaneously according to Eq. (8).

PDE Model: A Conservation Law for Locusts
We construct a continuous-time, mean-field model for the density of locusts. As outlined in Section 1.2, we write a continuous function of space and time R(x, t) for the density of available resources. Similarly, we write S(x, t) and M (x, t) for the density of stationary and moving locusts, respectively. See Table 1 for comparison with the variables of the agent-based model. These densities are governed by the partial differential equations which describe the feeding, switching, and movement behaviors on a macroscopic scale. The rate of decrease of R is proportional to the density of stationary locusts and available resources as described in Section 1.2. The constant of proportionality is given by the foraging rate λ. As in the ABM, locust foraging efficiency decreases as resources decrease. Note that the food R is decreasing in time at each spatial point x. The rate of change of S is determined wholly by the switching behavior. Here, the decrease of S represents the switching of locusts from stationary to moving with a rate dependent on R through k sm (R). Similarly, S increases as locusts switch from moving to stationary with rate k ms (R). See equation (2) for the functional forms of k sm , k ms . The same terms with opposite signs contribute to changes in M . The term vM x in the equation for M represents the marching of moving locusts to the right with the individual speed v. This spatial derivative makes the third equation into a standard transport equation.
A full list of all parameters appears in Table 2.
We consider initial conditions with resources that are a positive constant R + for large x; that is, R(x, 0) has lim x→∞ R(x, 0) = R + . We assume initial locust densities S(x, 0), M (x, 0) are smooth (continuous with continuous derivative). For biologically reasonable choices of such initial conditions, all solutions are guaranteed to remain positive, continuous, and finite by standard PDE theory of semi-linear conservation laws [28].
Finally, since the switching terms are of opposite signs in the S and M equations, we have mathematically guaranteed a conservation law. In particular, the total number of locusts in our 1-meter cross section N = Numerical Simulations For direct numerical simulations of the PDE, we use a 4th-order Runge-Kutta method for the temporal derivative with step dt. By choosing dx = v · dt we approximate the spatial derivative by a simple shift of the discritized M (x) on the spatial grid. For additional accuracy, we implement these schemes using a split-step method, as in [29] for instance. All simulations of the PDE used Matlab.

Parameter Identification
We identify a range of values for biological parameters from a variety of sources including research papers, Australian government guides and reports (particularly the Australian Plague Locust Commission), and agricultural organizations. A list of input parameters and ranges can be found in Table 2. A list of observable outcomes can be found in Table 3.

Input Parameters
Empirical Estimates We estimate five parameters directly from empirical observations: the total number of locusts N in the cross section, the initial resource density R + , the speed v of an individual locust, and the two switching rates when no resources are present k sm (R = 0), k ms (0). We provide ranges for these parameters in the first five rows of Table 2.
The total number of locusts N in our model is the number of locusts in a 1-meter cross section as shown in Fig 1. We rely on Buhl et al [9] to estimate N . In [9, Fig 1], the authors present three profiles of locust density computed by counting locusts in frames of video of a marching locust band taken during field experiments. The authors fit exponential curves to these data, see [9, Fig 2], which yield exponential rates of decay of density in time. We use these rates to estimate the area under the density profiles by integrating a corresponding exponential function. This provides three estimates for the total number of locusts who passed under the camera, which range from 9300 to 15000.
Rather than a precise measurement, we consider this an estimate and acknowledge that it may be improved by more direct analysis of the underlying data in [9]. We believe it does capture the correct order of magnitude and so include only a modestly larger range in our table.
Typical resource densities R + come from Meat and Livestock Australia [30]. This resource indicates that pasture with vegetation between 4 − 10cm high is desirable for livestock grazing. It also converts this range to a vegetation density measured in units of kilograms green Dry Matter per hectare. (Note that this measure discounts the mass of water in the vegetation, sometime up to 80%. While locusts typically feed on live non-dry vegetation, its water content does not provide energy or nutrients. As a result our variable R reflects not the harvestable greenery but instead represents the locust-edible resources.) We convert units and arrive at the range given in our table.
We obtained the speed v of an individual marching locust from experimental measurements in [31] and unpublished field data [22,32]. The experiments were conducted with the desert locust, S. Gregaria, and results in a range of 0.0339 − 0.0532 m/sec. This range contains the estimate from field data for the Australian plague locust from Buhl [22,32]. The latter sources also provides a second (higher) estimate that accounts for hopping, a common behavior of the Australian plague locust. Buhl's observations also show an increase in an individual's speed (averaged over crawling and hopping) with increasing temperature. Our range in Table 2 spans all three of these estimates. Most other recorded observations of speed represent collective informationthe speed of the aggregate band -which we discuss in Section 1.5.2 below.
Constants α and β represent the proportion of locusts that switch from stationary to moving (and vice versa) on bare ground, R ≈ 0. One laboratory study [21] with S. Gregaria provides data from which we draw out a single estimate for β as follows. The authors record the probability of these transitions in a laboratory area with no food present. They construct probability distributions (depending on time) for these transitions and fit curves to these distributions, see [21,Fig 1]. They find an exponential best fit for the probability that a locust transitions from moving to stationary. The exponential rate represents a reasonable value for β, so we gather that β ≈ 0.368 sec -1 . We use this estimate to set a minimum value of 0.01 < β and provide an upper limit below. The same source does not provide an estimate for α because the authors find that the probability distribution for stationary to moving transitions is best described by a power law.
Instead, for α we rely on the unpublished field data of [32]. A similar procedure as above yields an estimate of α ≈ 0.56 sec -1 . We use this to set a maximal value of 1 > α and provide a lower limit below.
Additional Parameters The parameters below the horizontal line in Table 2  chaos of the swarm. We discuss each and conduct a parameter sensitivity analysis in Section 2.4.
Constants η and θ represent the proportion of locusts that begin/restart or stop marching in a resource-rich environment, R ≈ R + . To empirically measure these would require a detailed examination of locusts marching in natural plant cover. We are not aware of a situation where such a study of marching has been conducted in a setting with abundant food.
To choose a range for η we rely on our biological assumption that a locust is more likely to begin moving when there are fewer resources nearby; that is, η < α. (In our sensitivity analysis of Section 2.4, this results in the bound η/α < 1.) This assumption provides a lower bound for α and an upper bound for η. We choose 0 as a lower bound for η, since it seems conceivable that a hungry locust might be satisfied to remain near food indefinitely. The converse biological assumption, that a locust is less likely to stop moving when there are fewer resources nearby, leads us to conclude that β < θ. (In our sensitivity analysis of, this results in the bound 1 < θ/β.) This provides our upper limit for β and our lower bound for θ. We choose our upper limit for θ to be significantly larger than η, the comparable transition rate with nearby food. This encodes an assumption that the attraction of nearby food is stronger than its absence. Note that these bounds are contained in the conditions we listed after introducing k sm , k ms in Eq (2). Namely, these choices force the transition rates to be decreasing and increasing respectively.
The parameters γ and δ determine how sharply the transition rates k sm (R) and k ms (R) depend on resources R. Specifically, they are the rate of exponential decrease and increase, respectively. One of our primary claims is that γ and δ must be positive, otherwise the transition rates k sm and k ms would be constant. More specifically, one may deduce that γ, δ should be of the same magnitude as 1/R + , since the functions k sm (R) and k ms (R) are defined on the interval [0, R + ]. Using our range of R + values above, we obtain the ranges appearing in the table for γ and δ.
The individual foraging rate λ is difficult to estimate for two reasons. First, it represents an instantaneous rate of change while most data on locust consumption is averaged over days or weeks, as in [34]. We found finer measurements of feeding in [33], where rates are averaged over ten-minute intervals. After unit conversions, we estimated a range of consumption rates on the order of 10 −8 − 10 −6 grams/(locust·sec). However, these rates are measured in a laboratory setting where locusts are provided with abundant resources to feed. This highlights a second difficulty in estimating λ; the data do not account for search times and so may represents a "feeding rate" rather than a foraging rate. Other factors such as digestion times complicate the matter further. With such persistent uncertainty, we allow a large range for λ and explore it thoroughly in our sensitivity analysis of Section 2.4.
Example Values Throughout the remainder of the text we illustrate our results using the set of example parameter values appearing in the second column from the right in Table 2. These values produce in both models a density profile consistent with observed locust bands. We selected these values using insight gleaned from our parameter sensitivity analysis, for details see the end of Section 2.4.

Collective Observables -Model Outcomes
We consider five measures of collective behavior. Table 3 provides an empirical range for each, estimated in the following paragraphs from data in the literature.
We approximated the collective speed c of the band from observations in [4,10,32]. Authors of [10] observed that bands moved between 36 − 92 meters per day (in "green grass"). [10, Table 4] estimates the times of day during which marching was observed, with a range of 3 to 7 total hours per day. We computed averages over these time intervals and converted units to obtain a range. In Clark [4], bands of locusts were observed for periods of an hour during daily marching. Both [4,32] report ranges of average band speeds overlapping with the range we computed above. Our Table 2 shows the union of all three ranges with rounding. Measuring this observable in our models is straightforward. In simulations of either the ABM or PDE we compute the mean position (or center of mass) of the locust band. Tracking the speed of the center of mass gives us the mean speed of the band. Additionally, analysis of the PDE model yields an explicit formula of for c with no need for simulations, details in Section 2.2.
The density of locust-edible food resources left behind by a band R − does not appear to be well studied. Wright [34] makes a careful study of leftover grain fit for human consumption; however, data are reported after threshing and processing and does not describe the amount of remaining green matter edible by locusts. An alternative approach to understanding R − could be to use [30] which suggests that a low range of green dry matter in pastures is 40 − 100 grams per square meter. This low range of green dry matter inhibits vegetation regrowth, increases erosion hazards, and is insufficient for grazing livestock. We emphasize that there is no data suggesting that a marching locust band leaves a field with leftover vegetation in this range. In particular, this provides us with an upper range only since some of the vegetation left behind may be inedible, even for voracious locusts. Thus we arrive at a lower bound of zero for R − . To measure the resources left behind in our models, we take a spatial average over the part of our domain to the left of the band of locusts.
The maximum locust density P = max (S + M ) in a band is taken from [10, Table  1]. We used the range of estimates observed for III and IV instars. This range is in line with the data of [34] who estimated a maximum density of 4000 locusts per square meter. In [9], the authors observed maximum densities ranging from 600 − 1200 locusts per square meter. We expect that these densities lie in (and just out of) the lower end of our range because the studies of [9] were conducted on bare ground with no vegetation while, typically, locusts aggregate into denser bands in lush vegetation, as observed in [4] for instance. The maximum density of a band in our models is measured simply by adding the components S and M and taking the maximal value.
The width W of the band, measured along the direction of motion, is taken primarily from Hunter et al [10]. Hunter et al measured the widths of bands by walking from the front into the band until "marching was no longer seen". Estimates from other sources fall in line with part of the range found in Hunter. For instance, 30 − 140 meters in [34] or 50 − 200 meters in [9]. We attribute the large range of band widths in [10] to the fact that these observations come from bands with a variety of sizes, as can also be seen by the large range for maximum densities in the same data set.
Measuring band width W in our model is not entirely straightforward as we cannot simply observe where "marching [is] no longer seen", as in [10]. Marching refers to a consistent movement of locusts with a preferred direction determined by alignment with their nearby neighbors. Since our models assume that locusts are always highly aligned, we rely on the locust density to determine where marching occurs. Experimental data and modeling work in [35] suggest that locusts in a group with a density greater than 20 locusts per square meter are likely to be highly aligned. We thus take W to be the length of the spatial interval where our density profile measures above the threshold of 20 locusts/m 2 .
This threshold definition of width W is biological and observable but it is not a good quantitative measure of the shape of a density distribution. For instance, consider a distribution with a maximum density less than the threshold density. This distribution will always measure W = 0 regardless of if it is very wide with a large total mass or if it is narrow with a much smaller mass. In other words, W does not scale with the total number of locusts in our band. We therefore introduce a second notion of width for use in comparing the shapes of bands with different total masses. A natural choice is the standard deviation of locust positions. We denote our standard deviation width by W σ and use it particularly in our sensitivity analysis of Section 2.4. Unfortunately, there is no general correspondence between our two notions of width W and W σ . Even for a fixed mass, one can construct distributions with different shapes and broad ranges of W σ while keeping W constant. For a given parameter set and varying mass we do compute W and make some a posteriori comparisons below.
The skewness Σ of a distribution is the third central moment and measures the distribution's symmetry about its mean. When Σ = 0 the distribution is symmetric while Σ > 0 suggests the distribution is leaning to the right with a longer tail on the left. (We acknowledge that this is the opposite of the standard convention.) Any exponential distribution e −Ax has skewness Σ = 2. Since [9] has demonstrated that an exponential fits well the locust density behind the peak, we consider 2 as a physically realistic upper bound. Including the sharp increase and maximum density at the front of the band will decrease skewness suggesting that we might expect values in the range 1 < Σ < 2.
Collective Observables for Example Values The example parameter values produce rather realistic collective outcomes with collective speed c = 0.0053 m/sec, remaining edible resources R − = 0.002 g/m, maximum density P = 1296 locusts/m 2 , threshold width W = 18.6 m (and standard deviation width W σ = 3.68 m), and skewness Σ = 1.78. Each of these values lies nearly within the ranges of Table 3. A small exception is the threshold width W , which is less than twelve units outside a wide range of several hundred units. Secondarily, we remind the reader of our difficulty in estimating the remaining resources. We interpret the small value R − = 002 g/m 2 to mean that in our models bands of locusts eat essentially all of the edible vegetable matter. We do not claim that they leave behind no vegetation at all.

ABM -Numerical Results
Typical behavior for the agent-based model is a transient period followed by a traveling pulse shape. During the transient period, the locust histogram variables S n,m and M n,m evolve to an equilibrium profile that moves with constant speed, each with stochastic variation at each time step. The duration of the transient period and shape of the equilibrium profile vary depending on biological input parameters, while the level of stochastic noise depends primarily on the size of ∆t. We explored a refinement of ∆t from 1 sec to 0.1 sec and observed similar behavior with decreasing levels of noise. In all results presented in this section we use ∆t = 1 sec and our example values from Table 2 for all biological parameters. The shape of the traveling pulse may be seen in Figure 4. The final histogram of locusts per spatial grid point at time t = 15000 appears in Fig 4a. A time-averaged pulse shape appears in Fig 4b. We construct this smooth density profile by averaging histograms for all time steps after the end of transients, in this case approximately t = 7500. Both plots show corresponding resource levels. The resources left behind R − after the pulse has completely passed depends primarily on the foraging rate constant λ. Shape of the traveling pulse profile also depends on λ but also on a complex combination of parameters in the stationary-moving transition probabilities p sm , p ms . For more detail on how the model depends upon parameters, see Section 2.4.
Qualitative and quantitative observations suggest that the tail of the density distribution of a hopper band is roughly exponential in shape [4,9,10]. Results from our agent-based model largely agree. We fit an exponential curve e a+bx to the tail of our average traveling pulse and obtained a Note the initial transients appearing as curves near t = 0, after which all each path appears piece-wise linear with either positive or negative slope corresponding to when the given locust was in a moving or stationary state. Each locust spends some time ahead of the mean and some time behind it, reminiscent of the "leap-frogging" behavior noted in [4].
convert the independent variable in the exponential from space in our numerical data to time in the empirical data. Since the pulse travels with with constant speed c, we have x = ct and our converted exponential is e a+bct with bc = 1.50 × 10 −3 , compared with exponential rates on the order of 10 −2 in [9].)

PDE -Hopper bands as traveling waves
Hopper bands require R-dependent switching To demonstrate the importance of the R-dependence in the switching rates k sm , k ms , we first consider a simplification of our model. Suppose that these switching rates are constant (k sm ≡ α, k ms ≡ β). We mathematically determine the long-time behavior of solutions to this simplified problem in S1 Appendix. For any locust density solution ρ = S + M , the center of mass moves to the right with a speed that approaches v α α+β as t → ∞. This is consistent with our search for traveling-wave solutions. However, we also find that the asymptotic standard deviation W σ ∼ √ t so that solutions spread diffusively for all time. In other words, no Existence of traveling wave solutions Returning to the main case with R-dependent switching, we show existence and development of hopper bands as traveling wave solutions to the PDE (9). A traveling wave is a solution with a fixed profile that propagates to the right with a constant speed c. S2 Appendix includes a mathematical analysis of traveling wave solutions. Numerical simulations suggest that these traveling waves organize all long-time dynamics of the model. That is, all solutions with our initial condition appear to converge to a traveling wave. Biologically, we conclude that a typical initial distribution of locusts aggregates into a coherent hopper band. The solid blue curves in Figure 5a show snapshots of the asymmetrical traveling wave created by R-dependent switching rates. Figures 5b and 5c compare the width and maximum density of the profiles for Locust density profiles with R-dependent (solid blue line) and R-independent (dashed gray line) switching rates. Each profile evolves from the same initial condition. Figure 5a shows snapshots of the density profiles over distance and time for both types of switching rates. Figure 5b illustrates the width of the bands where color represents a locust density greater than 20. Figure 5c displays the peak density of each pulse over time.
switching rates with and without resource dependence. In Figure 5b, colored regions correspond to a locust density greater than 20 locusts/m with gray and blue corresponding to R-independent and R-dependent switching rates, respectively. As the locust band without R-dependent switching progresses, the width of the gray region increases in time, showing diffusive spreading. On the other hand, the width of the locust band with R-dependent switching (blue) remains constant over time.
Additionally, the locust band with R-dependent switching reaches a constant height as seen in Fig 5c (blue). In contrast, the maximum locust density with R-independent switching rates decreases over time as locusts spread out (dashed gray).
Selection of principles By viewing hopper bands as traveling waves, our existence proof also determines a relationship between the total number (or total mass) N of locusts in our 1-meter cross section, the average band speed c, and the initial and remaining resources R + and R − . In S3 Appendix we show that these four variables must satisfy an explicit equation (32), (33) for any traveling wave. One consequence is that our model exhibits a selection mechanism whereby the average band velocity and the remaining resources are determined by the number of locusts in the band and the initial resource level. These explicit equations are illustrated in Fig 6. Each subfigure shows curves on which R + is constant (level curves). Plotting these in the N, c-plane (mass vs. speed), we obtain Fig 6a. (Here each curve is parameterized by R − .) Note that the curves appear monotone: speed c increases as a function of mass N . Biologically, this is what one expects; a larger swarm consumes food more quickly and moves on at a faster average pace. In Fig 6b, we plot the same level curves in the R − , c-plane (remaining resources vs. speed). (Now each curve is parameterized by N .) Again the curves are monotone but we now see that speed c decreases as a function of remaining resources R − . Here we also observe that the speed is much more sensitive when the remaining resources are very small. In S3 Appendix, we use the explicit formulas (32), (33) to prove the monotonicity of speed as a function of input parameters.

Agreement between ABM and PDE
We evaluate agreement between our two models by comparing the collective observables of Table 3. We divide these into two groups: the shape of the band as characterized by maximum density P, width W , and skewness Σ; and the mean speed c and remaining resources R − , which we consider the more agriculturally relevant.

ABM simulations and PDE analysis
The quantities c and R − can be determined for the PDE model via the traveling wave analysis of the last section. This analysis results in explicit formulas (32), (33). Substituting input parameters total mass N and initial resources R + , one can calculate exact results for c and R − . These relationships are represented by level curves in Fig 6, for details see Section 2.2. We ran direct numerical simulations of the ABM for selected values of the total mass N = 5000, 5250, 5500, 6000, 7500, 10000, 12500, 15000. In each simulation we used our example values for all other biological parameters. We ran each simulation for 2.5 × 10 4 time steps with ∆t = 1 for a final end time of 25000 sec and confirmed that the simulation reached the end of transients. We measured the collective speed c and remaining resources behind the band R − for each simulation. The resulting values agree with the explicit formulas to within 1% and are shown in Figure 6 (orange dots).
Direct simulation of both models We used direct numerical simulation of both models to evaluate their agreement on the basis of the shape characteristics maximum density P, standard deviation width W σ , and skewness Σ.
We ran both models for nt = 2 × 10 5 time steps using our example parameters and a range for the foraging rate λ so that −8 < log(λ) < −4. For each value of λ, we plot the shape characteristics in Fig 7. For the PDE, we measure the shape characteristics of the final output density profile. For the ABM, we measure the shape characteristics of a time-averaged density profile (as constructed in Fig 4b). The plots in Fig 7 are the result of continuation in the parameter log(λ). We begin with log(λ) = −4 and chose initial conditions computed from independent simulations of each model. For each value of log(λ) the algorithm proceeds as follows: We run both models for nt time steps; measure P, W σ , and Σ; choose new spatial grids for each model based on the value of W σ ; increase log(λ) by 0.1; and use the current output as the next initial condition. Practically speaking, the interval −8 < log(λ) < −4 is in fact covered by three such continuations originating at −4 and −7. Note that our numerical scheme begins to reach its limits as log(λ) approaches −8 because there the evolution of the profile shape is so slow that it requires very long computation times to reach equilibrium. This is also why we do not cover the full range of log(λ) explored in the next section.
To visually examine the agreement in the profiles, see Finally, these profiles also provide insight into the possible shapes of density profiles Comparison of the peak, width, and skewness of profiles from the PDE (blue line) and the ABM (orange circle), both obtained through direct numerical simulation for 2 × 10 5 time steps. Each shape observable is measured from the final numerical output for the PDE and from an time-averaged output the ABM. Longer simulations with 10 6 time steps, for ABM (gold x) and PDE (gold dot) show little evolution in the profile for longer times. and away from our example value of log(λ) = −5. Immediately, we notice that the remaining resources R − behind the pulse decrease quickly as foraging rate λ grows, confirming intuition. Next, the shape also varies dramatically as can be seen by noting that the horizontal axes in each row have vastly different scales. In particular, the profiles in the top row are short and wide while the middle row is narrow and tall, all having the same total number of locusts. The bottom row reveals a transition where the resources are nearly all depleted behind the pulse, leading to wide asymmetrical profile as observed in the field [10]. We carry out a more rigorous study of how the model responds to changes in the input parameters in the next section.

Parameter Sensitivity Analysis
The sensitivity of the model to its parameters was examined by computing Sobol indices [36] for several biologically observable quantities (see Table 3) with samples from the parameter space chosen via Saltelli's extension of the Sobol sequence [37,38]. Sobol indices represent a global, variance-based sensitivity analysis for nonlinear models that  Table 2 for parameter definition and ranges; this analysis was run using 4,400,000 samples from the given ranges. All log functions are base-10. First-order indices neglect all interactions with other parameters while total-order indices measure sensitivity through all higher-order interactions. Max 95% confidence intervals for the response variables was 0.01 for the first-order indices, 0.049 for total-order. has become extremely popular in recent years for examining the performance of mathematical models in the context of data (e.g., [39,40]). One of its strengths is the ability to calculate not just first-order (one-at-a-time) parameter sensitivity, but also second-order (two-at-a-time) and total-order (all possible combinations of parameters that include the given parameter) indices [38]. All indices are normalized by the variance of the output variable. Here, we will focus on the first-order and total-order indices, and note that the presence of higher-order interactions between the parameters can be inferred by comparing differences between these two. Scalar output quantities for our model (our collective observables) were all chosen with respect to the asymptotic traveling wave solution of the PDE model and are calculated by solving analytically for this solution. The observables chosen are the speed of the traveling wave c, the density of remaining resources R − , the peak (maximum) density of the wave profile, the width of the profile measured by its standard deviation W σ , and the skewness Σ of the profile. Section 1.5.2 provides physically relevant ranges for these observables from empirical studies.
In the case of switching parameters α, β, θ, and η, sensitivity to the ratios θ/β and η/α and the ratio difference ∆ = α/β − η/θ were used rather than the parameters themselves. One reason for this choice is to guarantee existence of a traveling wave solution; existence is guaranteed whenever ∆ > 0. Note that we also would like η/α < 1 and θ/β > 1 so that k sm is a decreasing function of resources R and k ms is an increasing function of resources. Additionally, these two conditions imply that ∆ > 0 so there is consistency between these constraints. Another reason for using these ratios lies in mathematical interpretation: ∆ is a measure of the difference in asymptotic switching rates behind the pulse (small R, α/β) and ahead of the pulse (large R, η/θ), and the two other ratios η/α (θ/β) describe how much the stationary to moving (moving to stationary) switching rates depend on R. More specifically, as these ratios approach 1 the switching rate changes little as R increases, while η/α close to zero or θ/β large implies a relatively large change in the switching rate. With these ratios and a value for β (chosen because we have some biological data for β), all four parameters in the ratios are uniquely determined.
Results are shown in Fig 9. All bars are stacked with each color corresponding to a different observable; reading across the parameters, the length of like colors can be compared. Critically, the parameter sensitivity is with respect to the range of parameter values given in the table included with Fig 9. These ranges were chosen to represent both biologically expected values (when information about these values could be obtained) and the necessary conditions for a traveling wave solution.
One immediate observation concerning the Sobol sensitivity analysis in Fig 9 is that log 10 (λ) and log 10 (∆) have a large effect on the density of resources asymptotically left behind the locust band as a fraction of the starting density (R − /R + ) and the ratio of the traveling wave velocity to the marching speed of a locust (c/v) respectively. Max density, pulse width as measured by standard deviation, and skewness also depend heavily on these two parameters. This is in fact unsurprising since these variables have by far the largest sample space range in terms of order of magnitude, and for this reason are the only ones examined on a log scale while all other parameters are on a linear scale. To explain this discrepancy, we remind the reader that our chosen sampling  Table 2 and represent 5% of all the points sampled for the Sobol analysis, chosen randomly. The red dot on the λ axis represents the example parameter set described in Table 2. ranges represent our uncertainty about the value that the parameters should take on in nature given all the information we were able to find in the biological literature. Our conclusion with this analysis then is that the model is in fact sensitive to this level of uncertainty in log 10 (λ) and log 10 (∆), and we should seek to narrow down the possibilities given what we know about observable, biological characteristics of the traveling locust band generated by our parameter choices (Table 3). Through the following numerical analysis of our sample data, we do just that.
To begin, we further illustrate the effect of varying log 10 (λ) and log 10 (∆) on the fraction of resources remaining R − /R + (in Fig 10) and on the ratio of the average pulse speed compared to the speed of a moving locust c/v (in Fig 11). In each figure, we plot a uniform random subset of the sample points used in the Sobol sensitivity analysis for the purpose of down-sampling the image and better visualizing sparse regions in the parameter space -it is qualitatively the same when using all sample points from the Sobol analysis.
Inspecting Fig 10 we note that, generally, at small λ a majority of resources persist after the locust front has passed while at large λ, the majority of resources are consumed. The red dot on the λ axis represents the example parameter set described in Table 2 which we believe to be a relatively feasible choice of parameter values in the context of the biological data about the observables in Table 3. We acknowledge that this appears to suggest the locusts leave behind no vegetation at all, but remember that our variable R represents locust-edible resources -there may be dry plant matter left behind that even locusts would not consume.  Table 2 and represent 5% of all the points sampled for the Sobol analysis, chosen randomly. The red dot represents the example parameter set described in Table 2.
Locust swarms observed in the field have a characteristic sharp rise at the beginning of the front and an exponential decay in the tail, see [9] for a quantitative analysis. This observation suggests that the skewness Σ is positive and less than or approximately equal to 2 (see Table 3). Figure 12 investigates the relationship between skewness Σ, foraging rate λ, and the difference of ratios ∆. For λ < 10 −7 , most values of Σ are negative, indicating an unrealistic density profile leaning to the left. As λ increases from 10 −7 to 10 −4 , Σ increases and clusters around 2. A smattering of points appear with Σ > 2 but these all correspond to profiles with unbiologically large maximum locust densities, as demonstrated by Fig 12b which only shows profiles with maximum locust density <10,000.
To identify a set of parameter inputs that would produce a density profile with observable quantities matching those found in the literature (see Section 1.5.2), we finally sorted the data underlying these figures and conditioned on desirable observable properties as specified in Table 3. This resulted in the example parameters specified in Table 2, with context provided by Figs 10,11,and 12. The results of the model run with these parameters can be seen in the figures included within the previous results sections.

Discussion
We present two minimal models for hopper bands of the Australian plague locust and demonstrate that resource consumption can mediate pulse formation. In these models all locusts are aligned and are either stationary (and feeding) or moving. Our agent-based model (ABM) tracks the locations, state, and resource consumption of individuals. In tandem, our partial differential equation (PDE) model represents the mean-field of the ABM. Both models generically form pulses as long as the transition rate from stationary to moving states is diminished by the presence of resources and/or the transition from moving to stationary states is enhanced by the presence of resources.
The agent-based model demonstrates that a locust's position within a travelling pulse circulates randomly. The PDE model provides a theoretical framework for proving the existence of traveling pulses. Additionally, this framework enables us to conduct an in-depth sensitivity analysis of the resulting pulse to the input parameters. Both models are capable of producing pulses whose shape and speed are consistent with field observations.
The ABM and the PDE each allow us to examine different facets of the problem. The ABM is easy to simulate and directly relates to microscopic field observations. It  Table 2 and represent 5% (in the case of 12a) and 50% (in the case of 12b) of all the points sampled for the Sobol analysis, chosen randomly. The red dot represents the example parameter set described in Table 2.
allows us to capture pulse formation, reproduces the stochastic variation seen in the field, and lets us track individual locusts, which perform random walks within the band. The PDE facilitates analysis of macroscopic behavior of the band including mean speed, total resource consumption, maximum locust density, and pulse width and skewness. Additionally, in the absence of resource dependence in the stationary-moving transition rates, we show that coherent pulses do not form -instead initial conditions spread indefinitely. These models are consistent in the following sense: the characteristics of pulses in the ABM, when averaged over many realizations, correlate precisely with the densities in the PDE model. We are fortunate that there is a healthy literature addressing the behavior of the Australian plague locust and the shape and speed of observed bands [4,9,10,14,32,33]. We have used these studies to estimate ranges for the microscopic parameters in our models. Some of these parameters (such as individual marching speed) have been carefully measured yielding narrow ranges. Others (notably the individual foraging rate) can only be deduced to lie within a range of several decades. Using these ranges, we analyze the sensitivity of a pulse's characteristics to changes in these input parameters. Sampling parameter values from within these ranges, we examine the resulting speed, remaining resources, and pulse peak, width, and skewness of over 4.4 × 10 6 traveling pulse profiles. Sobol sensitivity analysis quantifies the change in these characteristics as a function of the change in each input parameter. Guided by this analysis, we are able to identify a set of parameters that produce pulses consistent with those observed in the field. We conclude that our model is capable of reproducing field observations and that resource-dependent transitions are a consistent explanation for the formation and geometry of traveling pulses in locust hopper bands.
The testbed we have developed for modeling pulse formation in locust hopper bands has several natural extensions. First, we could incorporate variations and motion of locusts transverse to the primary direction of propagation. This two-dimensional model might capture the curving of pulses often seen in the field. Second, we could incorporate social interactions of locusts (notably alignment, attraction and collision avoidance). Lastly variations in resource density could be included. These extensions could help explain the variety of morphologies observed in hopper bands of Australian plague locusts and other species.
4 Supplemental Material S1 Appendix Resource-Independence: The Telegrapher's Equation. One of our primary assertions is that, in order for the locust population to form a coherent traveling pulse, the two transition rates k sm , k ms must depend on the resource density R. To demonstrate this we will estimate asymptotically the large-time behavior for the population densities in the case with constant transition rates.
Suppose k sm ≡ α and k ms ≡ β (which is equivalent to setting η = α and θ = β), then the equations governing the population densities become These equations are linear with constant coefficients and can be solved by a variety of means. Physically, they correspond to the probability densities associated to an agent-based model of random switching between stationary and moving states. In this interpretation, α is the probability that an agent transitions from stationary to moving and β is the probability of transition from moving to stationary. As such, we can identify this model as a variant of the Telegrapher's Equation.
Following previous work (see [41] and references therein), we use the method of moments to determine the large time behavior. We take initial conditions corresponding to a starting at the origin, where S 0 (M 0 ) is the probability is initially stationary (moving) and δ is the Dirac δ-function. We choose S 0 + M 0 = 1, reflecting the fact that this is a probability density. On average, an agent spends α α+β of its time moving. This leads us to conclude that its average velocity is c ≡ v α α+β . This motivates a change of variables ξ = x − ct which yields In this co-moving frame an agent now moves left with speed c (in the S state) and right with speed c − v (in the M state) but is never stationary. Solutions to this PDE correspond to probability distributions which can be characterized by their moments. We define the nth moment Multiplying (12) by ξ n and integrating yields the equations A similar calculation yield new initial conditions S n (0) = S 0 δ n,0 , M n (0) = M 0 δ n,0 where δ p,q is the Kronecker δ-function. When n = 0, we obtain the equations governing S 0 and M 0 the probability of a being in state S (or M ) at time t, The solution is As t → ∞, S 0 (t) ∼ β α+β and M 0 (t) ∼ α α+β as previously advertised. The first moments S 1 , M 1 , when divided by S 0 and M 0 , are the centers of mass for the corresponding the probability distributions S(ξ, t), M (ξ, t). To find these, we solve dS 1 dt = −αS 1 + βM 1 − cS 0 (t) which is guaranteed to be less than zero for some c as long as dksm dR ≤ 0 and dkms dR ≥ 0 and they are not both zero.
The ρ-nullcline is a vertical line occurring at R = R * where R * satisfies K(R * ) = 0. We must ensure that there is an interval of R + values such that R * ∈ (0, R + ). Note that K(0) > 0, K(∞) < 0, and K (R) < 0. By the continuity of P as R → ∞, we can apply the Intermediate Value Theorem and guarantee a unique R * ∈ (0, R + ) for any large enough choice of R + .
Fixing R + sufficiently large, we proceed with an invariant region argument, see Figure 13. Define the rectangle A = {0 < R ≤ R * , 0 < ρ < ρ * } for an arbitrary ρ * to be determined below. Note that region A is invariant as ξ decreases. This is simply due to the fact that R ξ > 0 and ρ ξ ≥ 0 on all of A. Therefore, any trajectory intersecting the ρ-nullcline {R = R * } must remain in region A as ξ decreases. By the Poincaré-Bendixson Theorem, the trajectory must terminate on some point (R − , 0) as ξ → −∞.
A B Fig 13. The (R, ρ) phase plane with semi-invariant regions A, B bounded by dotted lines including the ρ-nullcline (gold), sample arrows for the vector field (red), and a heteroclinic from (R − , 0) to (R + , 0) (blue).
We define region B by choosing the upper boundary as a line segment parallel to and ε-above the stable eigenspace of the linearization of (30) at (R + , 0) for some small ε > 0. That is, (Note that we may now choose the upper bound of region A to be ρ * = (R * ).) All that remains is to show that the stable manifold W s of (R + , 0) is contained within region B until it exits (as ξ decreases) through the vertical ρ-nullcline. Since W s is locally tangent to the stable eigenspace near (R + , 0), we know that it is contained in region B as ξ → ∞. Also W s cannot end, as ξ → −∞, on any of the equilibria composing the lower boundary {ρ = 0} of B. This is simply due to the fact that ρ ξ < 0 on the interior of region B so we must have that ρ increases as ξ decreases. Since R ξ > 0 on the vertical line segment between (R + , 0) and (R + , ), it is impossible for W s to exit B there as ξ decreases. Finally, W s cannot exit along because the slope of the vector field ρ ξ R ξ = 1 R c 1 − c K(R) at any point (R, ρ) is greater than the slope of . This follows from the fact that K (R) < 0 and that R < R + on all of . In fact, we have shown that, as ξ decreases, no trajectory may exit region B through any part of its boundary other than the vertical ρ-nullcline.

S3
Appendix Formulas for N, c, R + , R − . Following the existence result in S2 Appendix, we obtain explicit formulas that relate N, c, R + , and R − . Given any two of these variables and model parameter values, the following equations determine the other two variables: where We prove that equivalent equations hold for the nondimensionalized model (26).
Theorem 2. Given a traveling wave solution to (26). Then the quantities N, c, R + , R − satisfy where I 1 , I 2 are given by (34) with the nondimensional versions of k sm , k ms .

Proof.
A traveling wave solution satisfies the ODE (30). Since R ξ > 0 in the first quadrant of the phase plane, we can write ρ as a function of R along any heteroclinic. Thus dρ dξ = dρ dR dR dξ so integrating along a heteroclinic, we have We also have where K(R) is given in (31). Therefore we have proved equation (35).
Dividing the equation for R ξ in (30) by R and integrating the left hand side, we have Meanwhile, the right hand side gives us In fact, these equalities can be used to prove monotonicity of the mass-speed relation.
First, we show ds dR − < 0. Taking the derivative of (37), we get We will show the numerator, call it I, is negative. Dividing I by k sm (R − ) · k ms (R − ), we get The integrand ksm(R) ksm(R − ) − kms(R) kms(R − ) is less than zero for R − < R < R + . To see this, note that at R = R − this integrand is 0. Also we know that k sm is decreasing in R and k ms is increasing in R, so the first term decreases and the second term (without the negative sign) increases. Then the integrand is indeed negative for all R > R − .
Second, we will show ds dN > 0. Differentiating (37) and (38) respectively, we get and dS dN = dR − dN Setting these equal to each other, we obtain Substituting (43) into either (41) or (42) We have already shown that I < 0. Because R + > R − , the numerator is negative. The denominator is also negative, so have shown that ds dN > 0.