Agent-based and continuous models of hopper bands for the Australian plague locust: How resource consumption mediates pulse formation and geometry
Andrew J. Bernoff, Michael Culshaw-Maurer, Rebecca A. Everett, Maryann E. Hohn, W. Christopher Strickland, Jasper Weinburd
AAgent-based and continuous models of hopper bands for theAustralian plague locust: How resource consumptionmediates pulse formation and geometry
Andrew J. Bernoff , Michael Culshaw-Maurer , Rebecca A. Everett , Maryann E.Hohn , W. Christopher Strickland , Jasper Weinburd Department of Mathematics, Harvey Mudd College, Claremont, CA, USA Departments of Entomology and Nematology/Evolution and Ecology, University ofCalifornia Davis, Davis, CA, USA Department of Mathematics and Statistics, Haverford College, Haverford, PA, USA Mathematics Department, Pomona College, Claremont, CA, USA Department of Mathematics and Department of Ecology & Evolutionary Biology,University of Tennessee, Knoxville, TN, USAAll authors contributed equally to this work.* [email protected] 26, 2020 1/37 a r X i v : . [ n li n . AO ] F e b bstract Locusts are significant agricultural pests. Under favorable environmental conditionsflightless juveniles may aggregate into coherent, aligned swarms referred to as hopperbands. These bands are often observed as a propagating wave having a dense front withrapidly decreasing density in the wake. A tantalizing and common observation is thatthese fronts slow and steepen in the presence of green vegetation. This suggests thecollective motion of the band is mediated by resource consumption. Our goal is tomodel and quantify this effect. We focus on the Australian plague locust, for whichexcellent field and experimental data is available. Exploiting the alignment of locusts inhopper 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 ofindividuals and a partial differential equation model that describes locust density. Inboth these models, locust are either stationary (and feeding) or moving. Resourcesdecrease with feeding. The rate at which locusts transition between moving andstationary (and vice versa) is enhanced (diminished) by resource abundance. This effectproves essential to the formation, shape, and speed of locust hopper bands in ourmodels. From the biological literature we estimate ranges for the ten input parametersof our models. Sobol sensitivity analysis yields insight into how the band’s collectivecharacteristics vary with changes in the input parameters. By examining 4.4 millionparameter combinations, we identify biologically consistent parameters that reproducefield observations. We thus demonstrate that resource-dependent behavior can explainthe density distribution observed in locust hopper bands. This work suggests thatfeeding behaviors should be an intrinsic part of future modeling efforts.
Author summary
Locusts aggregate in swarms that threaten agriculture worldwide. Initially theseaggregations form as aligned groups, known as hopper bands, whose individualsalternate between marching and paused (associated with feeding) states. TheAustralian plague locust (for which there are excellent field studies) forms widecrescent-shaped bands with a high density at the front where locusts slow in uneatenvegetation. The density of locusts rapidly decreases behind the front where the majorityof food has been consumed.Most models of collective behavior focus on social interactions as the key organizingprinciple. We demonstrate that the formation of locust bands may be driven byresource consumption. Our first model treats each locust as an individual agent withprobabilistic rules governing motion and feeding. Our second model describes locustdensity with deterministic differential equations. We use biological observations ofindividual behavior and collective band shape to identify numerical values for the modelparameters and conduct a sensitivity analysis of outcomes to parameter changes. Ourmodels are capable of reproducing the characteristics observed in the field. Moreover,they provide insight into how resource availability influences collective locust behaviorthat may eventually aid in disrupting the formation of locust bands, mitigatingagricultural losses.
Introduction
Locusts are a significant agricultural pest in parts of Africa, Asia, Central and SouthAmerica, and Australia. They aggregate in large groups with as many as billions ofindividuals that move collectively, consuming large quantities of vegetation [1, 2].February 26, 2020 2/37ollective movement occurs in both nymphal and adult stages of development and isassociated with an epigenetic phase change from a solitary to a gregarious social statewhich is mediated by conspecific density and abiotic factors [1, 3–6]. Flightless nymphsmarch along the ground in aligned groups, often through agricultural systems wherethey cause significant crop damage as they feed and advance [4, 7, 8]. Some species, suchas the brown locust
Locustana pardalina , form intertwining streams of relativelyhomogeneous density [1, 2, 8]. By contrast the Australian plague locusts
Chortoicetesterminifera form wide, crescent-shaped bands that contain a high density in front and arapidly decreasing density behind [4, 9, 10]. Clark [4] notes:
The structure of bands varies according to the type of pasture through whichthey are passing. In areas of low cover containing plenty of green feed, bandsdevelop well-marked fronts in which the majority of hoppers may beconcentrated. In areas lacking green feed, bands lose their dense fronts andextend to form long streams, frequently exhibiting marked differences indensity throughout.
As bands of
C. terminifera move through a field of low pasture, they create a sharptransition from undamaged vegetation in front of the band to significant defoliationimmediately behind the band, see a schematic in Fig 1 or aerial photographs [10, Fig2], [11, Fig 1], [12, Fig 9], and [13]. In natural systems,
C. terminifera tend to consumeone of several species of grasses; in agricultural systems, they tend to eat primarilypasture and sometimes early stage winter cereals [14].
Fig 1.
Schematic of a traveling pulse of locusts. The Australian plague locust formsbroad hopper bands that propagate through vegetation in the direction perpendicular tothe aggregate structure [4, 9, 10]. The cross-sectional density profile is a traveling pulse ,with a steep leading edge (right) and shallower decay behind (left) that is roughlyexponentially decreasing in density [9]. Aerial photographs [10, Fig 2] show a notablecontrast between the verdant green of the unperturbed crops in front of the band andthe lifeless brown in the pulse’s wake. The one meter wide strip above represents thedimensions we use to model locust movement in a single dimension, as described below.The Australian plague locust
C. terminifera is the most common locust species onthe Australian continent. For ease, we henceforth refer to
C. terminifera simply as“locust”. Outbreaks of locust nymphs emerge as the result of a pattern of rainfall,vegetation growth, and drought [11, 15] which promotes breeding, hatching, crowding,and gregarization [3]. Gregarious nymphs form hopper bands of aligned individuals,which march distances from tens to hundreds of meters in a single day [4]. LocustsFebruary 26, 2020 3/37roceed through five nymphal stages, called instars , with marching behavior beginningduring the second instar [4]. Throughout these phases of life, hoppers consume largequantities of green biomass with an individual eating one third to one half of its bodyweight per day [14]. Approximately four weeks after eggs hatch, locusts reach adulthoodand are then capable of forming even more destructive and highly mobile flyingswarms [16].This study focuses on collective marching in hopper bands, which dominates thebehavior of gregarious locust nymphs in the third and fourth instars. Temperature andsunlight dictate a daily cycle of behavior with basking in the morning, roosting atmidday, and active periods of collective marching and feeding for up to nine hours whentemperatures are in an optimal range ( ∼ ◦ C) [10]. During these periods of collectivemarching, individuals crawl and hop across the ground in nearly the same direction astheir neighbors, due to social interactions [17, 18]. When individuals at the front of aband encounter available food resources, they stop and feed (see [4, 10] for qualitativeobservations and [19] for quantitative experimental results). Immediately after feeding,locusts exhibit a post-prandial quiescent period whose duration increases with theamount consumed [19–21]. Locusts farther back in the band may continue to moveforward, eventually passing those that stopped. This creates a “leap-frogging” typemotion with a cycling of individuals in the dense front of the band. Clark [4] describesthis behavior:
Those hoppers behind the front were in places which had been partly orwholly eaten out, and thus lacked the same stimulus of food to stop them. Astheir average rate of progress was greater than that of the hoppers in thefront, they tended to overtake them, becoming in turn slowed down in theirprogress by the presence of food.
Thus, individual motion during marching depends on individuals stopping to feed andconsequently on local resource density. We hypothesize that this effect mediates thecoherence and persistence of hopper bands with a dense front [4, 10] as well as thecharacteristic cross-sectional density distribution documented in [9].To test our hypothesis we conduct an in-depth modeling study concentrating on theinteraction of pause-and-go motion with food resources. We assume that hoppers marchin an aligned band through a field of finite resources, which is depleted as the locustsstop to feed. We develop a model for the probability of movement or stopping as afunction of resource availability. We construct and analyze in tandem an agent-basedmodel (ABM), which tracks individual locusts, and a partial differential equation (PDE)model, which considers mean-field densities. Both models produce traveling pulse typesolutions that are consistent with the detailed field observations of Buhl et al [9]. TheABM is easily simulated, allows us to track individuals within the swarm, and capturesthe natural stochasticity of a biological process. In contrast, the PDE produces smoothsolutions and lends itself to analysis and a detailed characterization of how observableoutcomes, such as mean band speed, cross-sectional density profile, and density ofresources left unconsumed in the wake, are related to the model’s parameters.Previous modeling efforts have considered both agent-based and continuous models,see [1] for an excellent overview of locust models. The majority of these have focused onsocial behavior – notably alignment, attraction, and repulsion with respect toconspecifics [18, 22–27]. Many of the agent-based models consider the pause-and-gobehavior of locusts [18, 26, 27], and other insects [28]. Continuous models have been usedto study transitions between stationary and moving states [29, 30] and gregarization [31].Foraging has been modeled in an agent-based framework [32] and resource distributioneffects on peak density has been posed as an energy minimization problem [33]. Othercontinuous models explicitly include food resources having animal movement depend ona combination of aggregation and gradient sensing (chemotaxis in many, startingFebruary 26, 2020 4/37ith [34], or “herbivory-taxis” in [35, 36], for instance). These studies find that travelinganimal bands are the result of a balance between attraction to food and inter-animaldispersal, bearing some qualitative resemblance to the results presented here. However,locusts in the present model do not sense resource gradients (instead, direction isprescribed implicitly by social alignment) and the corresponding mathematicalequations are distinct from the well-studied equations of chemotaxis.We are aware of no models of locust band movement that incorporate foragingbehavior or food resources. Previous studies such as [23, 27] suggest that the formationof sharp asymmetric fronts may be explained solely by social forces. By contrast, ourmain conclusion is that foraging and resource-mediated stationary/moving transitionsproduce pulse-shaped density profiles, supporting the observations of hopper bands withdense fronts and the inferences on foraging of Clark [4] and Hunter et al. [10]. A furtherstrength of our model is that it quantitatively reproduces the observed density profilesof [9] from biologically realistic parameters.In Section 1 we construct our two models beginning with biological and simplifyingassumptions, and ending with parameter identification from empirical field data inTable 2. Section 2 contains our results, namely that both models produce a travelingpulse in locust density precisely when the locusts’ stationary/moving transitions aredependent upon the amount of nearby resources. Evidence consists of numericalsimulations for the ABM, mathematical traveling wave analysis for the PDE, and arobust sensitivity analysis of the models to changes in the input parameters. In Section3 we revisit our main findings and outline extensions of this work incorporating morebiological complexity. Our appendices in the Supporting Information containmathematical analysis and proofs substantiating results for the PDE.
We outline our assumptions for the modeling framework. Our models are minimal inthe sense that we include only the effects necessary to investigate the main question:
Can resource-dependent locust behavior drive the formation of a dense front and thepropagation of hopper bands? • We assume that resources (food) can only decrease, since locusts feed much morequickly than vegetation grows. Moreover, resources are identical so that they canbe characterized by a single variable. Prior to locust arrival, we assume availableresources have a spatially uniform density. • We model only the part of the daily cycle dominated by collective movement.During a typical day, a hopper band has one or two periods of collectivemovement (marching) totaling up to nine hours. The remainder of the day isspent resting (basking and roosting) [1, 4, 8–10]. • We assume hopper bands consist of flightless nymphs that are behaviorallyidentical in all regards. Bands often include a mix of two instars (e.g. II and IIIinstars or III and IV instars) which behave qualitatively similarly with laterinstars being larger, eating more, and moving more quickly. • We assume individuals move parallel to one another, creating a constant directionof movement for the entire band. Locusts are known to align their direction ofmovement with their nearest neighbors and may align with environmental cuessuch as wind or the location of the sun [2, 4, 8].February 26, 2020 5/37
We model behavior in a narrow strip aligned to the direction of movement, asshown in Fig 1. For dimensional consistency of the model, we assume thetransverse width to be 1 meter. • We assume that each individual is either stationary or moving. Further, onlystationary locusts feed while moving locusts propagate forward with a constantspeed that represents an average of crawling and hopping. • We assume that locusts feed continuously when they are stationary. In fact,locusts eat a meal and then remain sedentary during a post-prandialperiod [19–21]. While biologically different, these processes are mathematicallyanalogous and we believe including such a delay in the model is unlikely tosignificantly alter our results.Furthermore, we make additional assumptions on the rate of transitions betweenmoving and stationary that are supported by empirical observation, although theycombine and simplify multiple locust behaviors. • We assume that locusts transition back and forth between stationary and movingstates at a rate depending solely on the resources nearby.Notably, we have not included any explicit social interaction between locusts;interaction is mediated solely through the consumption of resources. Socialinteraction plays a well-document role in the aggregation, alignment, andmarching of hopper bands, see [1] for instance. By modeling one spatial dimensiononly, we implicitly include the social tendency of locusts to align their direction ofmotion with neighbors as demonstrated in [17]. We do not focus on socialinteractions simply because our primary goal is to investigate the effect of linkingresource consumption with pause-and-go motion on hopper band morphologies. • We assume the transition rate from moving to stationary is positive and increasesas the resource density increases.Field observations [4, 10] and laboratory experiments [19] have shown thatindividuals stop marching to eat when they encounter resources in their path.While we assume resources have a uniform local density, the reality on the groundis that a locust is more likely to encounter an edible plant, and thus stop to feed,when the resource density is high. • We assume the transition rate from stationary to moving is positive anddecreasing with resource density.This behavior is consistent with foraging theories, such as the simple mechanismsillustrated in [37] where insects are likely to leave a patch of resources before thepoint of diminishing returns. The
Marginal Value Theorem [38] quantifies thisbehavior: if an energy cost assigned to foraging is proportional to resource density,then when local resource densities drop below a critical level it costs less energyper unit resource to move on in search of higher density resources. Additionally,there is a second, more subtle behavior behind this assumption. Locusts thatbecome stationary are assumed to have consumed resources. After feeding, locustsexhibit a post-prandial period of inactivity which extends in proportion with theamount consumed [19–21]. Our assumption about this transition rate reflects alonger period of inactivity when resources are plentiful and larger amounts aretherefore consumed.Foreshadowing our results, only one of these two transition rates must dependstrictly on local resource availability for our model to produce coherent travelingpulse-type density profiles akin to observed hopper bands with dense fronts.February 26, 2020 6/37
These transition processes are completely memoryless, which implies that locustsexperience neither hunger nor satiation.The biological reality is that feeding behavior is complex, see [39] for a review.Locust hunger has been well documented in other species [19, 40]. Since, in ourmodel, locusts are traveling through a field of relatively plentiful resources wesuggest that most locusts do not experience starvation (i.e. no sustenance for 24hours as in some experiments).We remind the reader that our goal in this study is to demonstrate thatresource-dependent behavior is sufficient for the formation and propagation of hopperbands with a coherent dense front. We acknowledge that the efficacy of this model maybe improved by adding social interactions – such as alignment, attraction, and repulsion.Additionally, we believe these additions, particularly that of alignment, would play apivotal role when modeling locust behavior in two-dimensions, as in [23, 27].
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. Table1 compares their independent and state variables and Table 2 lists their common modelparameters.In the ABM, space and time lie on a discrete, evenly spaced lattice ( x n , t m ) while inthe PDE space and time ( x, t ) are continuous. In both models, S and M denote thenumber or density of stationary and moving locusts respectively. For the ABM, thenumber 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 resourcedensity in the continuous model is R ( x, t ).We assume that the group rate of feeding is proportional to the product of thestationary 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 locustsconsume resources. This implies that an individual locust’s foraging efficiency decreasesas resources become scarcer at their location. This is not an explicit implementation ofthe Marginal Value Theorem but fits the general concept of foraging efficiency within apatch decreasing due to searching time, not satiation by the forager [38]: as theresources at a location are eaten, locusts have difficulty locating the next unit toconsume, reducing the overall rate of resource consumption at that location. We willrefer to λ as the foraging rate , as it reflects both feeding and foraging efficiency.We model the stationary-moving transitions as a Markov (memoryless) process. Forthe PDE model, this yields a rate at which the population of stationary locuststransitions to moving, and vice versa. This assumption ignores the transition historyand hunger (as discussed above) of any individual locust, which is justifiable on thetimescale of the collective motion (hours). We use exponentially saturating functions ofresources as illustrated in Fig 2. The stationary to moving rate is denoted k sm while themoving to stationary rate is called k ms . Specifically, k sm ( R ) = η − ( η − α ) e − γR , k ms ( R ) = θ − ( θ − β ) e − δR , (2)where γ, δ >
0, 0 ≤ β ≤ θ , and 0 < η ≤ α . The conditions on the parameters guaranteethat k sm ( R ) is a decreasing function and k ms ( R ) is a increasing function of R . (Most ofFebruary 26, 2020 7/37he analytical results concerning the PDE model hold for any choice of monotoneswitching rates - see S2 Appendix for details.) Fig 2.
Transition rates for stationary to moving k sm (gold) and moving to stationary k ms (purple) with α = 0 . β = 0 . γ = 0 . δ = 0 . η = 0 . θ = 0 . k sm , k ms appear as coefficients in growth anddecay terms in the differential equations. In the ABM we use a stochastic version ofthese transitions. At each time step, locusts switch from stationary to moving via atransition probability p sm and from moving to stationary via p ms , both of which arefunctions of R n,m . The smooth transition rates k ms and k sm can be understood to bederived from these probabilities as the time step ∆ t approaches zero. Assuming ∆ t issmall yields the following approximations, p sm ( R n,m ) ≈ k sm ( R n,m )∆ t, p ms ( R n,m ) ≈ k ms ( R n,m )∆ t. (3)This is equivalent to assuming that each locust undergoes only a single transition in anygiven time step. Biologically, these transition probabilities can be estimated fromintermittent motion observed in the laboratory [41] or the field [42]. These observationssuggest that transitions occur on a timescale of a few second. Additionally locusts alsoexhibit a post-prandial quiescence which may last several minutes, particularly after alarge meal [19–21]. These timescales are much shorter than the period of collectivemarching (hours) which justifies our original approximation that the process isMarkovian (memoryless). In using ranges for α and β , we aim to allow for naturalvariation between the two species for which there is data.February 26, 2020 8/37gent-Based Units Continuous Units DescriptionModel Model x n L x L position (along direction of motion) t m T t T time S n,m C S ( x, t ) P number/density of stationary locusts M n,m C M ( x, t ) P number/density of moving locusts R n,m Q R ( x, t ) Q edible resource density Table 1.
Independent and dependent variables appearing in the agent-based andpartial differential equation models. Units are L = length [meters], T = time [seconds], C = number of locusts, P = locust density [number/(meter) ], and Q = resourcedensity [grams/(meter) ]. We now describe the details and implementation of our agent-based model (ABM) whichencodes the behavior of each individual locust. The temporal evolution of the ABMmay be thought of as a probabilistic cellular automaton. The model is one dimensionalin 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), andthe spatial availability of resources (food). Locust position and the spatial distributionof 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 onestep on the lattice per 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 statevariable where σ i = 1 when the locust is moving and 0 otherwise. The motion of thelocusts can now be expressed succinctly as X i ( t m +1 ) = X i ( t m ) + σ i ( t m ) v ∆ t = X i ( t m ) + σ i ( t m )∆ x, (4)where we have applied the value of the state variable at t m throughout the interval oflength ∆ t . Note this artifice ensures that the values X i remain on the lattice for alltime t m .We model transitions between stationary and moving states with a discrete-timeMarkov process given via the probabilities in Eq (3). Thus, at time t m , each locust atposition x n has a probability p sm ( R n,m ) to switch from stationary σ i = 0 to moving σ i = 1 or a probability p ms ( R n,m ) to switch from moving to stationary.We define the histogram variables mentioned above by simply counting the numberof locusts in each state at each space-time grid point: M n,m = (cid:88) X n ( t m ) σ n ( t m ) = σ i = 1) locusts at ( x n , t m ) (5) S n,m = (cid:88) X n ( t m )(1 − σ n ( t m )) = σ i = 0) locusts at ( x n , t m ) . (6)We model the resources with a scalar variable R n,m which is defined as availablefood, 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 dR n,m dt = − λ S n,m ∆ x R n,m . (7)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 − λS n,m ∆ t ∆ x . (8)February 26, 2020 9/37iologically, this evolution implies that the resources in a patch of vegetation infestedby a group of stationary, feeding locusts will decrease by approximately half in anamount 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/ ( λ · R n, = R + , indicating an initially constant field of resources.Together with initial conditions, Eqs (3), (4), and (8) specify the evolutioncompletely. Our agent-based model then takes the form of three sequential, repeatingsteps 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, andresources in each location are also updated simultaneously according to Eq (8). We construct a continuous-time, mean-field model for the density of locusts. Asoutlined in Section 1.2, we write a continuous function of space and time R ( x, t ) for thedensity of available resources. Similarly, we write S ( x, t ) and M ( x, t ) for the density ofstationary and moving locusts, respectively. See Table 1 for comparison with thevariables of the agent-based model.These densities are governed by the partial differential equations R t = − λSRS t = − k sm S + k ms MM t = k sm S − k ms M − vM x x ∈ R , t ∈ [0 , ∞ ) , (9)which describe the feeding, switching, and movement behaviors on the scale of theaggregate band. The rate of decrease of R is proportional to the density of stationarylocusts and available resources as described in Section 1.2. The constant ofproportionality is given by the foraging rate λ . As in the ABM, locust foragingefficiency decreases as resources decrease. Note that the food R is decreasing in time ateach spatial point x . The rate of change of S is determined wholly by the switchingbehavior. Here, the decrease of S represents the switching of locusts from stationary tomoving with a rate dependent on R through k sm ( R ). Similarly, S increases as locustsswitch from moving to stationary with rate k ms ( R ). See Eq (2) for the functional formsof 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 withthe individual speed v . This spatial derivative makes the third equation into a standardtransport equation. A full list of all parameters appears in Table 2.We consider initial conditions with resources that are a positive constant R + forlarge x ; that is, R ( x,
0) has lim x →∞ R ( x,
0) = R + . We assume initial locust densities S ( x, , M ( x,
0) are non-negative and smooth (continuous with continuous derivative).For biologically reasonable choices of such initial conditions, all solutions are guaranteedto remain non-negative, continuous, and finite by standard quasilinear hyperbolicPDE [43].Finally, since the switching terms are of opposite signs in the S and M equations, wehave mathematically guaranteed a conservation law. In particular, the total number oflocusts in our 1-meter cross section N = (cid:82) ∞−∞ ( S + M ) dx is conserved.February 26, 2020 10/37 umerical Simulations For direct numerical simulations of the PDE, we use a4th-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 discretized M on the spatial grid. This is equivalent to a first-order upwind scheme because M ( x n ,t m +1 ) − M ( x n ,t m ) dt = − v M ( x n ,t m ) − M ( x n − ,t m ) dx = ⇒ M ( x n , t m +1 ) = M ( x n − , t m ) . For additional accuracy, we implement these schemes using a split-step method, asin [44] for instance. All simulations of the PDE used Matlab.
We identify a range of values for biological parameters from a variety of sourcesincluding research papers, Australian government guides and reports (particularly theAustralian Plague Locust Commission), and agricultural organizations. A list of inputparameters and ranges can be found in Table 2. A list of observable outcomes can befound in Table 3.
We estimate five parameters directly from empiricalobservations: the total number of locusts N in the cross section, the initial resourcedensity R + , the speed v of an individual locust, and the two switching rates when noresources are present k sm ( R = 0) , k ms (0). We provide ranges for these parameters inthe first five rows of Table 2.The total number of locusts N in our model is the number of locusts in a 1-metercross 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 inframes of video of a marching locust band taken during field experiments. The authorsfit exponential curves to these data, see [9, Fig 2], which yield exponential rates of decayof density in time. We use these rates to estimate the area under the density profiles byintegrating a corresponding exponential function. This provides three estimates for thetotal 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 thatit may be improved by more direct analysis of the underlying data in [9]. We believe itdoes capture the correct order of magnitude and so include only a modestly larger rangein our table.Typical resource densities R + come from Meat and Livestock Australia [45]. Thisresource indicates that pasture with vegetation between 4 − R reflects not the harvestable greenery but instead represents thelocust-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 experimentalmeasurements in [46] and unpublished field data [27, 42]. The experiments wereconducted with the desert locust, Schistocerca gregaria , and results in a range of0 . − . α and β represent the proportion of locusts that switch from stationary tomoving (and vice versa) on bare ground, R ≈
0. One laboratory study [26] with
S.gregaria provides data from which we draw out a single estimate for β as follows. Theauthors record the probability of these transitions in a laboratory area with no foodpresent. They construct probability distributions (depending on time) for thesetransitions and fit curves to these distributions, see [26, Fig 1]. They find anexponential best fit for the probability that a locust transitions from moving tostationary. The exponential rate represents a reasonable value for β , so we gather that β ≈ .
368 sec -1 . We use this estimate to set a minimum value of 0 . < β and providean upper limit below. The same source does not provide an estimate for α because theauthors find that the probability distribution for stationary to moving transitions is bestdescribed by a power law.Instead, for α we rely on the unpublished field data of [42] for C. terminifera . Asimilar procedure as above yields an estimate of α ≈ .
56 sec -1 . We use this to set amaximal value of 1 > α and provide a lower limit below. In using ranges for α and β ,we aim to allow for natural variation between the two species for which there is data. Description Units Min Max Example Source N total number locusts in strip C/L R + resource density in front of band Q
120 250 200 [45] v individual marching speed L/T α S → M transition rate for R = 0 1 /T η β M → S transition rate for R = 0 1 /T θ η S → M transition rate, large R /T α θ M → S transition rate, large R /T β γ exponent of S → M transition 1 /Q δ exponent of M → S transition 1 /Q λ individual foraging rate 1 /T P − − − [47] Table 2.
Estimates of biological parameters for both models. Parameters above thehorizontal line are estimated from empirical observations, with explanations in text.Parameters below the horizontal are estimated from collective information and modelbehavior. Units are L = length [meters], T = time [seconds], C = number of locusts, P = locust density [number/(meter) ], and Q = resource density [grams/(meter) ]. Additional Parameters
The parameters below the horizontal line in Table 2 do notall have readily available estimates in the literature; likely because the individualinformation encoded in these parameters is difficult to measure empirically amid thechaos of the swarm. We discuss each and conduct a parameter sensitivity analysis inSection 2.4.Constants η and θ represent the proportion of locusts that begin/restart or stopmarching in a resource-rich environment, R ≈ R + . To empirically measure these wouldrequire a detailed examination of locusts marching in natural plant cover. We are notaware of a situation where such a study of marching has been conducted in a settingwith abundant food.To choose a range for η we rely on our biological assumption that a locust is moreFebruary 26, 2020 12/37ikely to begin moving when there are fewer resources nearby; that is, η < α . (In oursensitivity analysis of Section 2.4, this results in the bound η/α < α and an upper bound for η . We choose 0 as a lower boundfor η, since it seems conceivable that a hungry locust might be satisfied to remain nearfood indefinitely. The converse biological assumption, that a locust is less likely to stopmoving when there are fewer resources nearby, leads us to conclude that β < θ. (In oursensitivity analysis of, this results in the bound 1 < θ/β .) This provides our upper limitfor β and our lower bound for θ . We choose our upper limit for θ to be significantlylarger than η , the comparable transition rate with nearby food. This encodes anassumption that the attraction of nearby food is stronger than its absence. Note thatthese bounds are contained in the conditions we listed after introducing k sm , k ms inEq (2). Namely, these choices force the transition rates to be decreasing and increasingrespectively.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 decreaseand 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, onemay 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 + valuesabove, we obtain the ranges appearing in the table for γ and δ. The individual foraging rate λ is difficult to estimate for two reasons. First, itrepresents an instantaneous rate of change while most data on locust consumption isaveraged over days or weeks, as in [48]. We found finer measurements of feeding in [47],where rates are averaged over ten-minute intervals. After unit conversions, we estimateda range of consumption rates on the order of 10 − − − grams/(locust · sec). However,these rates are measured in a laboratory setting where locusts are provided withabundant resources to feed. This highlights a second difficulty in estimating λ ; the labdata does not account for search times and so may represents a “consumption rate”rather than a foraging rate. To explain, recall that our ABM places a locust at a gridpoint which represents a rectangle of physical space with dimensions ∆ x × . Alocust may need to move within this small rectangle to find an individual plant suitablefor feeding. Since we track only the resource density in that local rectangle, this searchtime is simply accounted for in the foraging rate. Other factors such as digestion timesand the post-prandial rest period complicate the matter further. With such persistentuncertainty, we allow a large range for λ and explore it thoroughly in our sensitivityanalysis of Section 2.4. Example Values
Throughout the remainder of the text we illustrate our resultsusing the set of example parameter values appearing in the second column from theright in Table 2. These values produce in both models a density profile consistent withobserved locust bands. We selected these values using insight gleaned from ourparameter sensitivity analysis, for details see the end of Section 2.4.
We consider five measures of collective behavior. Table 3 provides an empirical range foreach, estimated in the following paragraphs from data in the literature.We approximated the collective speed c of the band from observations in [4, 10, 42].Authors of [10] observed that bands moved between 36 −
92 meters per day (in “greengrass”). [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 timeintervals and converted units to obtain a range. In Clark [4], bands of locusts wereobserved for periods of an hour during daily marching. Both [4, 42] report ranges ofFebruary 26, 2020 13/37verage band speeds overlapping with the range we computed above. Our Table 2 showsthe union of all three ranges with rounding. Measuring this observable in our models isstraightforward. In simulations of either the ABM or PDE we compute the meanposition (or center of mass) of the locust band. Tracking the speed of the center of massgives us the mean speed of the band. Additionally, analysis of the PDE model yields anexplicit 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 notappear to be well studied. Wright [48] makes a careful study of leftover grain fit forhuman consumption; however, data are reported after threshing and processing anddoes not describe the amount of remaining green matter edible by locusts. Analternative approach to understanding R − could be to use [45] which suggests that alow range of green dry matter in pastures is 40 −
100 grams per square meter. This lowrange of green dry matter inhibits vegetation regrowth, increases erosion hazards, and isinsufficient for grazing livestock. We emphasize that there is no data suggesting that amarching 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 maybe 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 thepart of our domain to the left of the band of locusts. Symbol Description Units Min Max Example CitationOutput c speed of collective band L/T . .
009 0.0053 [4, 10, 42] R − remaining resource density Q P
950 4280 1296 [9, 10, 48] W threshold width of profile L
30 500 18.6 [9, 10, 48]Σ skewness of locust profile 1 1 2 1.78 [9]
Table 3.
Collective observables with ranges based on field research. Units are L =length [meters], T = time [seconds], C = number of locusts, P = locust density[number/(meter) ], and Q = resource density [grams/(meter) ]. Note that skewness (Σ)is nondimensional.The maximum locust density P = max ( S + M ) in a band is taken from [10, Table1]. We used the range of estimates observed for III and IV instars. This range is in linewith the data of [48] who estimated a maximum density of 4000 locusts per squaremeter. In [9], the authors observed maximum densities ranging from 600 − S and M and taking the maximal value.The width W of the band, measured parallel to the direction of motion, is takenprimarily from Hunter et al [10]. Hunter et al measured the widths of bands by walkingfrom the front into the band until “marching was no longer seen”. Estimates from othersources fall in line with part of the range found in Hunter. For instance, 30 −
140 metersin [48] or 50 −
200 meters in [9]. We attribute the large range of band widths in [10] tothe fact that these observations come from bands with a variety of sizes, as can also beseen 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 cannotFebruary 26, 2020 14/37imply observe where “marching [is] no longer seen”, as in [10]. Marching refers to aconsistent movement of locusts with a preferred direction determined by alignment withtheir nearby neighbors. Since our models assume that locusts are always highly aligned,we rely on the locust density to determine where marching occurs. Experimental dataand modeling work in [17] suggest that locusts in a group with a density greater than 20locusts per square meter are likely to be highly aligned. We thus take W to be thelength of the spatial interval where our density profile measures above the threshold of20 locusts/m .This threshold definition of width W is biological and observable but it is not a goodquantitative measure of the shape of a density distribution. For instance, consider adistribution with a maximum density less than the threshold density. This distributionwill always measure W = 0 regardless of if it is very wide with a large total mass or if itis narrow with a much smaller mass. In other words, W does not scale with the totalnumber of locusts in our band. We therefore introduce a second notion of width for usein comparing the shapes of bands with different total masses. A natural choice is thestandard 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 isno general correspondence between our two notions of width W and W σ . Even for afixed 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 docompute W and make some a posteriori comparisons below.The skewness Σ of a distribution is the third central moment and measures thedistribution’s symmetry about its mean. When Σ = 0 the distribution is symmetricwhile Σ > e − Ax has skewness Σ = 2. Since [9] has demonstrated that an exponentialfits well the locust density behind the peak, we consider 2 as a physically realistic upperbound. Including the sharp increase and maximum density at the front of the band willdecrease skewness suggesting that we might expect values in the range 1 < Σ < Collective Observables for Example Values
The example parameter valuesproduce rather realistic collective outcomes; each of them is very nearly in the rangeobtained from the literature, see Table 3. A small exception is the threshold width W ,which is less than twelve units outside a large range of several hundred units.Secondarily, we remind the reader of our difficulty in estimating the remaining resources.We interpret the small value R − = 0 .
002 g/m to mean that in our models bands oflocusts eat essentially all of the edible vegetable matter. We do not claim that theyleave behind no vegetation at all. Typical behavior for the agent-based model is a transient period followed by a travelingpulse 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 withstochastic variation at each time step. The duration of the transient period and shapeof the equilibrium profile vary depending on biological input parameters, while the levelof 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 allresults presented in this section we use ∆ t = 1 sec and our example values from Table 2for all biological parameters.February 26, 2020 15/37ig 3A shows the instantaneous speed of the mean position of all locusts over thecourse of 10000 sec. After an initial increase, the speed stabilizes around an average c = 5 . × − m/sec with a standard deviation of 0 . × − m/sec. Individual locustsmove according to a biased random walk around the mean position, as illustrated by thepaths of five sample locusts in Fig 3B. Note the brief period of transients visible asarcing curves near t = 0, after which the distance to the mean is given by a piece-wiselinear function for each locust. Intervals with positive slope v − c correspond to periodswhere the individual was moving, while negative slope − c indicates periods where theindividual was stationary and the mean position marched on ahead. Fig 3. (A) Speed of the mean position of all locusts (center of mass of the swarm).Note the initial increase followed by a sustained period of variation around the average c = 5 . × − m/sec. The standard deviation around c is 0 . × − m/sec aftertransients. (B) Paths of five sample locusts, each shown in a different color. Note theinitial transients appearing as curves near t = 0, after which all each path appearspiece-wise linear with either positive or negative slope corresponding to when the givenlocust was in a moving or stationary state. Each locust spends some time ahead of themean and some time behind it, reminiscent of the “leap-frogging” behavior noted in [4].The shape of the traveling pulse may be seen in Fig 4. The final histogram of locustsper spatial grid point at time t = 15000 appears in Fig 4A. A time-averaged pulse shapeappears in Fig 4B. We construct this smooth density profile by averaging histograms forall time steps after the end of transients, in this case approximately t = 7500. Bothplots show corresponding resource levels. The resources left behind R − after the pulsehas completely passed depends primarily on the foraging rate constant λ . Shape of thetraveling pulse profile also depends on λ but also on a complex combination ofparameters in the stationary-moving transition probabilities p sm , p ms . For more detailon how the model depends upon parameters, see Section 2.4.Qualitative and quantitative observations suggest that the tail of the densitydistribution of a hopper band is roughly exponential in shape [4, 9, 10]. Results from ouragent-based model agree. We fit an exponential curve e a + bx to the tail of our averagetraveling pulse and obtained a = 4 . , b = 0 . ig 4. Output of the agent-based model with N = 7000 locust agents at time t = 15000 sec, re-centered so that the mean position of all locusts occurs at zero. (A)Shows the final state of the model including the number of locusts at each spatial gridpoint (orange) and the remaining resource density at each spatial grid point (dotted,green). Compares well with previously published data [9, Fig 1]. (B) Displays atime-average of model outputs taken after an amount of time to account for transients(in this case, approximately t = 7500). Gray shading indicates ± one standard deviationfrom the average locust density (blue) and resources (green). The tail of the pulseagrees well with an exponential least squares fit (gold).convert the independent variable in the exponential from space in our numerical data totime 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 . × − , compared withexponential rates on the order of 10 − in [9].) Hopper bands require R -dependent switching To demonstrate the importanceof the R -dependence in the switching rates k sm , k ms , we first consider a simplification ofour model. Suppose that these switching rates are constant ( k sm ≡ α , k ms ≡ β ). Wemathematically determine the long-time behavior of solutions to this simplified problemin S1 Appendix. For any locust density solution ρ = S + M , the center of mass movesto the right with a speed that approaches v αα + β as t → ∞ . This is consistent with oursearch for traveling-wave solutions. However, we also find that the asymptotic standarddeviation W σ ∼ √ t so that solutions spread diffusively for all time. In other words, noFebruary 26, 2020 17/37oherent hopper bands form in the long-time limit. Gray dashed lines in Fig 5A depictthis behavior, illustrating the decay of a locust density profile with resource-independentswitching rates. Fig 5.
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.(A) Shows snapshots of the density profiles over distance and time for both types ofswitching rates. (B) Illustrates the width of the bands where color represents a locustdensity greater than 20. (C) Displays the peak density of each pulse over time. Existence of traveling wave solutions
Returning to the main case with R -dependent switching, we show existence and development of hopper bands astraveling wave solutions to the PDE (9). A traveling wave is a solution with a fixedprofile that propagates right or left with a constant speed c . Since locusts move only tothe right in our model, we expect right-moving traveling waves and S2 Appendixincludes a mathematical analysis of these solutions. Numerical simulations suggest thatthese traveling waves organize all long-time dynamics of the model. That is, allsolutions with our initial condition appear to converge to a traveling wave. Biologically,we conclude that a typical initial distribution of locusts aggregates into a coherenthopper band. The solid blue curves in Fig 5A show snapshots of the asymmetricaltraveling wave created by R -dependent switching rates.Fig 5B and 5C compare the width and maximum density of the profiles for switchingrates with and without resource dependence. In Fig 5B, colored regions correspond to alocust density greater than 20 locusts/m with gray and blue corresponding to R -independent and R -dependent switching rates, respectively. As the locust bandwithout R -dependent switching progresses, the width of the gray region increases intime, 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 bandwith R -dependent switching reaches a constant height as seen in Fig 5C (blue). Incontrast, the maximum locust density with R -independent switching rates decreasesover time as locusts spread out (dashed gray).February 26, 2020 18/37 raveling waves dynamically select collective observables By viewing hopperbands as traveling waves, our existence proof also determines a relationship between thetotal number (or total mass) N of locusts in our 1-meter cross section, the average bandspeed c , and the initial and remaining resources R + and R − . In S3 Appendix we showthat these four variables must satisfy an explicit equation for any traveling wave. Oneconsequence is that our model exhibits a selection mechanism whereby the average bandvelocity and the remaining resources are determined by the number of locusts in theband and the initial resource level.These explicit equations are illustrated in Fig 6. Each subfigure shows curves onwhich 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 curvesappear monotone: speed c increases as a function of mass N . Biologically, this is whatone expects; a larger swarm consumes food more quickly and moves on at a fasteraverage pace. In Fig 6B, we plot the same level curves in the R − , c -plane (remainingresources vs. speed). (Now each curve is parameterized by N .) Again the curves aremonotone 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 remainingresources are very small. In S3 Appendix, we use the explicit formulas to prove themonotonicity of speed as a function of input parameters. Fig 6.
Level curves on which initial resources R + are constant (black curves) computedexplicitly from the analytic formulas of the PDE model and numerical data points(orange circle) generated by direct simulation of the ABM with N = 5000 , , , , , , , We evaluate agreement between our two models by comparing the collective observablesof Table 3. We divide these into two groups: the shape of the band as characterized bymaximum density P, width W , and skewness Σ; and the mean speed c and remainingresources R − , which we consider to be more agriculturally relevant. ABM simulations and PDE analysis
The quantities c and R − can be determinedfor the PDE model via the traveling wave analysis of the last section. This analysisresults in explicit formulas in S3 Appendix. Substituting input parameters total mass N and initial resources R + , one can calculate exact results for c and R − . Theserelationships 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 , , , , , , , . × time steps with ∆ t = 1 for a final end time of 25000 sec and confirmed that theFebruary 26, 2020 19/37imulation reached the end of transients. We measured the collective speed c andremaining resources behind the band R − for each simulation. The resulting values agreewith the explicit formulas to within 1% and are shown in Fig 6 (orange circle). Direct simulation of both models
We used direct numerical simulation of bothmodels to evaluate their agreement on the basis of the shape characteristics maximumdensity P, standard deviation width W σ , and skewness Σ. Fig 7.
Comparison of the peak, width, and skewness of profiles from the PDE (blueline) and the ABM (orange circle), both obtained through direct numerical simulationfor 2 × time steps. Each shape observable is measured from the final numericaloutput for the PDE and from a time-averaged output the ABM. Longer simulationswith 10 time steps, for ABM (gold x) and PDE (gold dot) show little evolution in theprofile for longer times. Note that the maximum density is higher for long simulationsof the ABM (gold x) because these represent a single instance, rather than an average.We ran both models for nt = 2 × time steps using our example parameters and arange for the foraging rate λ so that − < log( λ ) < −
4. For each value of λ , we plot theshape characteristics in Fig 7. For the PDE, we measure the shape characteristics of thefinal output density profile. For the ABM, we measure the shape characteristics of atime-averaged density profile (as constructed in Fig 4B). The plots in Fig 7 are theresult of continuation in the parameter log( λ ). We begin with 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 − < log( λ ) < − − −
7. Note that our numerical scheme begins toreach its limits as log( λ ) approaches − λ ) explored in the next section.To visually compare the profiles, see Fig 8. These six profiles are the result ofrunning each model for nt = 10 time steps with our example parameter values andselected log( λ ) = − . , − . , − .
2. First, note the strong agreement along each row.Second, a data point (gold dot and x) from each of these log( λ ) values is included in theplots of Fig 7. Since there is little difference between these data points and the rest inthe figure, which are the result of only 2 × time steps, we can conclude that theshape characteristics have reached near-equilibrium values. The gold x at log( λ ) = − . λ ) = −
5. Immediately, we notice that the remainingresources R − behind the pulse decrease quickly as foraging rate λ grows, confirmingintuition. Next, the shape also varies dramatically as can be seen by noting that thehorizontal axes in each row have vastly different scales. In particular, the profiles in thetop row are short and wide while the middle row is narrow and tall, all having the sametotal number of locusts. The bottom row reveals a transition where the resources arenearly all depleted behind the pulse, leading to wide asymmetrical profile as observed inthe field [10].We carry out a more rigorous study of how the model responds to changes in theinput parameters in the next section. The sensitivity of the model to its parameters was examined by computing Sobolindices [49] for several biologically observable quantities (see Table 3) with samples fromthe parameter space chosen via Saltelli’s extension of the Sobol sequence [50, 51]. Sobolindices represent a global, variance-based sensitivity analysis for nonlinear models thathas become extremely popular in recent years for examining the performance ofmathematical models in the context of data (e.g., [52, 53]). One of its strengths is theability to calculate not just first-order (one-at-a-time) parameter sensitivity, but alsosecond-order (two-at-a-time) and total-order (all possible combinations of parametersthat include the given parameter) indices [51]. All indices are normalized by thevariance of the output variable. Here, we will focus on the first-order and total-orderindices, and note that the presence of higher-order interactions between the parameterscan be inferred by comparing differences between these two.Scalar output quantities for our model (our collective observables) were all chosenwith respect to the asymptotic traveling wave solution of the PDE model and arecalculated by solving analytically for this solution. The observables chosen are thespeed 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 itsstandard deviation W σ , and the skewness Σ of the profile. Section 1.5.2 providesphysically 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 parametersthemselves. One reason for this choice is to guarantee existence of a traveling wavesolution; existence is guaranteed whenever ∆ >
0. Note that we also would like η/α < θ/β > k sm is a decreasing function of resources R and k ms is anincreasing function of resources. Additionally, these two conditions imply that ∆ > R , α/β ) and ahead of the pulse (large R , η/θ ),and the two other ratios η/α ( θ/β ) describe how much the stationary to movingFebruary 26, 2020 21/37moving to stationary) switching rates depend on R . More specifically, as these ratiosapproach 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 anda value for β (chosen because we have some biological data for β ), all four parameters inthe ratios are uniquely determined.Results are shown in Fig 9. All bars are stacked with each color corresponding to adifferent observable; reading across the parameters, the length of like colors can becompared. Critically, the parameter sensitivity is with respect to the range of parametervalues given in the table included with Fig 9. These ranges were chosen to representboth biologically expected values (when information about these values could beobtained) and the necessary conditions for a traveling wave solution.One immediate observation concerning the Sobol sensitivity analysis in Fig 9 is thatlog ( λ ) and log (∆) have a large effect on the collective observables of the pulse. Fig 8.
Model outputs from direct numerical simulation for 10 time steps. Densityprofiles from the PDE (blue, left) and histograms from the ABM (orange, right) forselected foraging rates log( λ ) = − . , − . , − .
2. For quick visual shape comparison, alloutputs shifted so that center of mass is x = 0. Each plot corresponds to a data point inFig 7 (gold x for ABM, gold dot for PDE). Note the differences in scale on thehorizontal axes in each plot.February 26, 2020 22/37 R + v β η / α θ / β l o g ( Δ ) γ δ l o g ( λ ) First-order indices max densitystd widthskewness N R + v β η / α θ / β l o g ( Δ ) γ δ l o g ( λ ) Total-order indices N R + v β η / α θ / β l g ( Δ Δ γ δ l g ( λ Δ c/vR + /R + N R + v β η / α θ / β l g ( Δ Δ γ δ l g ( λ Δ Value RangeN 5000-30000R + Fig 9.
Sensitivity of various traveling wave observables to model parameters (bars arestacked). See Table 2 for parameter definition and ranges; this analysis was run using4,400,000 samples from the given ranges. All log functions are base-10. First-orderindices neglect all interactions with other parameters while total-order indices measuresensitivity through all higher-order interactions. Max 95% confidence intervals for theresponse variables was 0.01 for the first-order indices, 0.049 for total-order.Recall that λ is the parameter encoding the foraging rate; ∆ is discussed in detailearlier in this section. The bottom row of Fig 9 shows the impact of these parameters Fig 10.
Scatter plot of (A) remaining resource fraction R − /R + and (B) fraction of thetraveling wave speed c over individual locust speed v as a function of the foraging rate λ and colored by ∆, the difference in asymptotic switching rates behind and ahead of thepulse. Points are taken from the parameter ranges in Table 2 and represent 5% of allthe points sampled for the Sobol analysis, chosen randomly. The red dots represent theexample parameter set described in Table 2.February 26, 2020 23/37n the density of resources asymptotically left behind the locust band as a fraction ofthe starting density ( R − /R + ) and the ratio of the traveling wave velocity to themarching speed of a locust ( c/v ) respectively. Max density, pulse width as measured bystandard deviation, and skewness also depend heavily on these two parameters as seenin the top row of Fig 9. This is in fact unsurprising since λ and ∆ have by far thelargest sample space range in terms of order of magnitude, and for this reason are theonly ones examined on a log scale while all other parameters are on a linear scale. Toexplain this discrepancy, we remind the reader that our chosen sampling rangesrepresent our uncertainty about the value that the parameters should take on in naturegiven all the information we were able to find in the biological literature. Ourconclusion with this analysis then is that the model is in fact sensitive to this level ofuncertainty in log ( λ ) and log (∆), and we should seek to narrow down thepossibilities given what we know about observable, biological characteristics of thetraveling locust band generated by our parameter choices (Table 3). Through thefollowing numerical analysis of our sample data, we do just that.To begin, we further illustrate the effect of varying log ( λ ) and log (∆) on thefraction of resources remaining R − /R + (in Fig 10A) and on the ratio of the averagepulse speed compared to the speed of a moving locust c/v (in Fig 10B). In each figure,we plot a uniform random subset of the sample points used in the Sobol sensitivityanalysis for the purpose of down-sampling the image and better visualizing sparseregions in the parameter space – it is qualitatively the same when using all samplepoints from the Sobol analysis.Inspecting Fig 10A we note that, generally, at small λ a majority of resources persistafter the locust front has passed while at large λ , the majority of resources areconsumed. The red dot on the λ axis represents the example parameter set described inTable 2 which we believe to be a relatively feasible choice of parameter values in thecontext of the biological data about the observables in Table 3. We acknowledge thatthis appears to suggest the locusts leave behind no vegetation at all, but remember thatour variable R represents locust-edible resources – there may be dry plant matter leftbehind that even locusts would not consume. Fig 11.
Skewness as a function of foraging rate ( λ ) and colored by ∆. Fig 11A isrepresentative of the entire sampled parameter space while Fig 11B shows only thepoints with peak wave amplitudes less than 10,000 locusts per square meter. Points aretaken from the parameter ranges in Table 2 and represent 5% (in the case of Fig 11A)and 50% (in the case of 11B) of all the points sampled for the Sobol analysis, chosenrandomly. 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 beginningof the front and an exponential decay in the tail, see [9] for a quantitative analysis. Thisobservation suggests that the skewness Σ is positive and less than or approximatelyFebruary 26, 2020 24/37qual to 2 (see Table 3). Fig 11 investigates the relationship between skewness Σ,foraging rate λ , and the difference of ratios ∆. For λ < − , most values of Σ arenegative, indicating an unrealistic density profile leaning to the left. As λ increases from10 − to 10 − , Σ increases and clusters around 2. A smattering of points appear withΣ > < We present two minimal models for hopper bands of the Australian plague locust anddemonstrate that resource consumption can mediate pulse formation. In these modelsall locusts are aligned and are either stationary (and feeding) or moving. Ouragent-based model (ABM) tracks the locations, state, and resource consumption ofindividuals. In tandem, our partial differential equation (PDE) model represents themean-field of the ABM. Both models generically form pulses as long as the transitionrate from stationary to moving states is diminished by the presence of resources and/orthe transition from moving to stationary states is enhanced by the presence of resources.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 observations at the scale ofindividual locusts. It captures pulse formation and propagation, reproduces thestochastic variation seen in the field, and lets us track individual locusts which performrandom walks within the band. The PDE model provides a theoretical framework forproving the existence of traveling pulses. This framework facilitates analysis of thecollective behavior of the band including mean speed, total resource consumption,maximum locust density, pulse width, and pulse skewness. In turn, this theory enablesus to conduct an in-depth sensitivity analysis of the pulse’s characteristics with respectto the input parameters. The two models are consistent in the following sense: thecharacteristics of pulses in the ABM, when averaged over many realizations, correlateprecisely with the densities in the PDE model.We are fortunate that there is a healthy literature addressing the behavior of theAustralian plague locust, notably the shape and speed of observed bands [4, 9, 10, 16, 47].We have used these studies to estimate ranges for the organism-level parameters in ourmodels. Some of these parameters (such as individual marching speed) have beencarefully measured yielding narrow ranges. Others (notably the individual foraging rate)can only be deduced to lie within a range of several orders of magnitude. Using thesebiologically plausible ranges, we analyze the sensitivity of a pulse’s characteristics tochanges in the input parameters. Sampling parameter values from within these ranges,we examine the resulting speed, remaining resources, and pulse peak, width, andskewness of over 4 . × traveling pulse profiles. Sobol sensitivity analysis quantifiesthe 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 produces pulsesconcordant with those observed in the field. We conclude that resource-dependenttransitions are a consistent explanation for the formation and geometry of travelingpulses in locust hopper bands.A reasonable question is whether a different mechanism can drive pulse formation orFebruary 26, 2020 25/37f the formation of pulses is enhanced by a combination of behaviors. Prior works, bothfor the Australian plague locust and for other locust species, investigate a variety ofsocial mechanisms for collective movement in hopper bands. In two agent-basedmodeling studies [23, 27], pulses are among a handful of aggregate band structuresobtained by varying the parameters that model individual locust behavior. Acontinuum approach in [30] finds traveling pulses in a PDE similar to our Eq (9) butwithout accounting for resources. Instead, social behavior is encoded via dependence onlocust density of both the transition rates and the speed. This is coarsely akin to ourmodel where resource-dependent transitions between moving and stationary states isnecessary for pulse formation, see S1 Appendix. However, a model with social behavioras the only driving factor does not account for the observations of Clark [4] andHunter [10] that C. terminifera manifests pulse-shaped bands with varying shapes basedon the surrounding vegetation. We believe that incorporating both social andresource-dependent behaviors will better reproduce field observations.Further experiments and field observations could help to elucidate the combinedroles of resources and social behavior in the formation of hopper bands. While there is aconsiderable literature on the social [3, 9] and feeding [15, 47] behavior of
C. terminifera ,less progress has been made in quantifying the effect of food on individuals in densebands exhibiting collective motion. One notable exception is the recent study of Dkhiliet al. [19]. Looking ahead, field data could be collected as video while a hopper bandmoved through lush vegetation, see the methods in [9]. With continuing advances inmotion tracking, for instance as employed in [28], one could collect time-series data oneach individual moving through the frame. From such data one could draw out theeffects of nearby vegetation, satiation or hunger, and local locust density onpause-and-go motion. In turn these processes could be modeled more thoroughly.We see our present models as a testbed upon which one may develop extensions thatcapture more of the complexity in locust hopper behavior. The most natural of theseextensions is to consider locusts’ social behavior, as discussed above. A second is toinclude stochastic, individual, and environmental variation. This could be incorporatedinto the agent-based model in order to examine the robustness of pulse characteristicswith respect to a distribution of individual marching speeds, or even large hops, asin [27]. For the PDE model, random variations in locust movement could naturally berepresented by a linear diffusion term. Thirdly, we could incorporate motion of locuststransverse to the primary direction of propagation. This two-dimensional model mightaim to capture the curving of the front of hopper bands often seen in the field. Lastly,large changes in resource density could be included to represent the band entering orexiting a lush field or pasture, with a view towards informing barrier control strategiesas discussed in [11, 19]. These extensions could help explain the variety of morphologiesand density profiles – including curving dense fronts, complex fingering, andlower-density columns – observed in hopper bands of the Australian plague locust andother species.
Acknowledgments
This collaboration began as part of a Mathematics Research Communities of theAmerican Mathematical Society on Agent-based Modeling in Biological and SocialSystems (2018), under National Science Foundation grant DMS-1321794. A follow-upvisit was also supported by a collaboration travel grant from the AMS MRC program.A subsequent research visit was supported by the Institute for Advanced Study SummerCollaborators Program. Jasper Weinburd gratefully acknowledges the support of anNSF Mathematical Sciences Postdoctoral Research Fellowship grant DMS-1902818.Christopher Strickland was supported in part by a grant from the Simons FoundationFebruary 26, 2020 26/37rant
References
1. Ariel G, Ayali A. Locust Collective Motion and Its Modeling. PLOSComputational Biology. 2015; p. 1–25. doi:10.1371/journal.pcbi.1004522.2. Uvarov B. Grasshoppers and Locusts. vol. 2. London, UK: Cambridge UniversityPress; 1977.3. Gray LJ, Sword GA, Anstey ML, Clissold FJ, Simpson SJ. Behavioural PhasePolyphenism in the Australian Plague Locust (
Chortoicetes Terminifera ).Biology Letters. 2009;5(3):306–309. doi:10.1098/rsbl.2008.0764.4. Clark LR. Behaviour of Swarm Hoppers of the Australian Plague locustChortoicetes terminifera (Walker). Commonwealth Scientific and IndustrialResearch Organization, Australia. 1949;245:5–26.5. Simpson SJ, McCaffery AR, H¨agele BF. A Behavioral Analysis of Phase Changein the Desert Locust. Biol Rev. 1999;74(4):461–480.6. Pener MP, Simpson SJ. Locust Phase Polyphenism: An Update. In: Advances inInsect Physiology. vol. 36. Elsevier; 2009. p. 1–272.7. Love G, Riwoe D. Economic costs and benefits of locust control in easternAustralia. Canberra: Australian Plague Locust Commission; 2005. 05.14.8. Symmons PM, Cressman K. Desert Locust Guidelines. Rome: UNFAO; 2001.9. Buhl J, Sword GA, Clissold FJ, Simpson SJ, Buhl J, Sword GA, et al. Groupstructure in locust migratory bands. Behavioral Ecology and Sociobiology.2011;65(2):265–273.10. Hunter DM, Mcculloch L, Spurgin PA. Aerial detection of nymphal bands of theAustralian plague locust (Chortoicetes terminifera (Walker)) (Orthoptera :Acrididae). Crop Protection. 2008;27:118–123. doi:10.1016/j.cropro.2007.04.016.11. Hunter DM. Advances in the Control of Locusts (Orthoptera: Acrididae) inEastern Australia: From Crop Protection to Preventive Control. AustralianJournal of Entomology. 2004;43(3):293–303. doi:10.1111/j.1326-6756.2004.00433.x.12. Silliman BR, McCoy MW, Angelini C, Holt RD, Griffin JN, van de Koppel J.Consumer Fronts, Global Change, and Runaway Collapse in Ecosystems. AnnualReview of Ecology, Evolution, and Systematics. 2013;44(1):503–538.doi:10.1146/annurev-ecolsys-110512-135753.13. Spratt W. Risk management of a major agricultural pest in Australia - plaguelocusts. The Australian Journal of Emargency Management. 2004;19(3):20–25.February 26, 2020 27/374. Walton CS, Hardwick L, Hanson J. Locusts in Queensland: Pest status reviewseries - Land protection. Queensland Government. 1994;78(February 2003):1193.doi:10.1094/PD-78-1193.15. Clissold FJ, Sanson GD, Read J. The paradoxical effects of nutrient ratios andsupply rates on an outbreaking insect herbivore, the Australian plague locust.Journal of Animal Ecology. 2006;75(4):1000–1013.doi:10.1111/j.1365-2656.2006.01122.x.16. Wright DE. Analysis of the development of major plagues of the Australianplague locust Chortoicetes terminifera (Walker) using a simulation model.Australian Journal of Ecology. 1987;12(4):423–437.doi:10.1111/j.1442-9993.1987.tb00959.x.17. Buhl J, Sumpter DJT, Couzin ID, Hale JJ, Despland E, Miller ER, et al. FromDisorder to Order in Marching Locusts. Science. 2006;312:1402–1406.doi:10.1126/science.1125142.18. Buhl J, Sword GA, Simpson SJ. Using field data to test locust migratory bandcollective movement models. Interface Focus. 2012;2:757–763.19. Dkhili J, Maeno KO, Idrissi Hassani LM, Ghaout S, Piou C. Effects of starvationand Vegetation Distribution on Locust Collective Motion. Journal of InsectBehavior. 2019;32(3):207–217. doi:10.1007/s10905-019-09727-8.20. Kennedy JS, Moorhouse JE. Laboratory observations on locust responses towind-borne grass odour. Entomologia Experimentalis et Applicata.1969;12(5):487–503. doi:10.1111/j.1570-7458.1969.tb02547.x.21. Moorhouse JE. Experimental analysis of the locomotor behaviour of Schistocercagregaria induced by odour. Journal of Insect Physiology. 1971;17(5):913–920.doi:10.1016/0022-1910(71)90107-7.22. Vicsek T, Czirok A, Ben Jacob E, Cohen I, Shochet O. Novel Type ofPhase-Transition in a System of Self-Driven Particles. Phys Rev Lett.1995;75(6):1226–1229.23. Dkhili J, Berger U, Idrissi Hassani LM, Ghaout S, Peters R, Piou C.Self-Organized Spatial Structures of Locust Groups Emerging from LocalInteraction. Ecological Modelling. 2017;361:26–40.doi:10.1016/j.ecolmodel.2017.07.020.24. Gr´egoire G, Chat´e H, Tu Y. Moving and Staying Together Without a Leader.Physica D. 2003;181(3-4):157–170.25. Romanczuk P, Couzin ID, Schimansky-Geier L. Collective Motion due toIndividual Escape and Pursuit Response. Physical Review Letters.2009;010602:1–4. doi:10.1103/PhysRevLett.102.010602.26. Ariel G, Ophir Y, Levi S, Ben-Jacob E, Ayali A. Individual pause-and-go motionis instrumental to the formation and maintenance of swarms of marching locustnymphs. PLOS One. 2014;9(7):e101636.27. Bach A. Exploring locust hopper bands emergent patterns using parallelcomputing. Universit´e Paul Sabatier Toulouse III; 2018.February 26, 2020 28/378. Nilsen C, Paige J, Warner O, Mayhew B, Sutley R, Lam M, et al. SocialAggregation in Pea Aphids: Experiment and Random Walk Modeling. PLOSOne. 2013;8(12):e83343.29. Navoret L. A two-species hydrodynamic model of particles interacting throughself-alignment. Mathematical Models and Methods in Applied Sciences.2013;23(06):1067–1098.30. Edelstein-Keshet L, Watmough J, Grunbaum D. Do travelling band solutionsdescribe cohesive swarms? An investigation for migratory locusts. Journal ofMathematical Biology. 1998;36:515–549. doi:10.1007/s002850050112.31. Topaz CM, D’Orsogna MR, Edelstein-Keshet L, Bernoff AJ. Locust Dynamics:Behavioral Phase Change and Swarming. PLOS Comp Bio. 2012;8(8):e1002642.32. Nonaka E, Holme P. Agent-Based Model Approach to Optimal Foraging inHeterogeneous Landscapes: Effects of Patch Clumpiness. Ecography.2007;30(6):777–788. doi:10.1111/j.2007.0906-7590.05148.x.33. Bernoff AJ, Topaz CM. Biological Aggregation Driven by Social andEnvironmental Factors: A Nonlocal Model and its Degenerate Cahn-HilliardApproximation. SIAM J Appl Dyn Sys. 2016;15(3):1528–1562.34. Keller EF, Segel LA. Traveling bands of chemotactic bacteria: A theoreticalanalysis. Journal of Theoretical Biology. 1971;30(2):235–248.doi:10.1016/0022-5193(71)90051-8.35. Gueron S, Liron N. A model of herd grazing as a travelling wave, chemotaxis andstability. Journal of Mathematical Biology. 1989;27(5):595–608.doi:10.1007/bf00288436.36. Lewis MA. Spatial Coupling of Plant and Herbivore Dynamics: The Contributionof Herbivore Dispersal to Transient and Persistent “Waves” of Damage.Theoretical Population Biology. 1994;45(3):277–312. doi:10.1006/tpbi.1994.1014.37. Bell WJ. Searching Behavior Patterns in Insects. Annual Review of Entomology.1990;35(1):447–467. doi:10.1146/annurev.en.35.010190.002311.38. Charnov EL. Optimal Foraging, the Marginal Value Theorem. TheoreticalPopulation Biology. 1976;9(2):129–136. doi:10.1016/0040-5809(76)90040-X.39. Simpson SJ, Raubenheimer D. The Hungry Locust. Advances in the Study ofBehavior. 2000; p. 1–44. doi:10.1016/s0065-3454(08)60102-3.40. Bernays EA, Chapman RF. Meal size in nymphs of
Locusta migratoria .Entomologia Experimentalis et Applicata. 1972;15(4):399–410.doi:10.1111/j.1570-7458.1972.tb00227.x.41. Bazazi S, Bartumeus F, Hale JJ, Couzin ID. Intermittent Motion in DesertLocusts: Behavioural Complexity in Simple Environments. PLoS ComputationalBiology. 2012;8(5):e1002498. doi:10.1371/journal.pcbi.1002498.42. Buhl J. Personal Communication; 2019.43. Evans LC. Partial differential equations. Providence, R.I.: AmericanMathematical Society; 2010.44. Bernoff AJ. Slowly varying fully nonlinear wavetrains in the Ginzburg-Landauequation. Physica D: Nonlinear Phenomena. 1988;30(3):363–381.February 26, 2020 29/375. Meat & Livestock Australia. Improving pasture use with the MLA Pasture Ruler;2014. Available from: https://mbfp.mla.com.au/pasture-growth/ .46. Bazazi S, Buhl J, Hale JJ, Anstey ML, Sword GA, Simpson SJ, et al. CollectiveMotion and Cannibalism in Locust Migratory Bands. Current Biology.2008;18(10):735–739. doi:10.1016/j.cub.2008.04.035.47. Bernays EA, Chapman RF. The role of food plants in the survival anddevelopment of Chortocietes terminifera (Walker) under drought conditions.Australian Journal of Zoology. 1973;21:575–592.48. Wright DE. Damage and loss in yield of wheat crops caused by the Australianplague locust,
Chortoicetes terminifera (Walker). Australian Journal ofExperimental Agriculture. 1986;26(5):613–618. doi:10.1071/EA9860613.49. Sobol IM. Global sensitivity indices for nonlinear mathematical models and theirMonte Carlo estimates. Math and Comput in Simul. 2001;55(1-3):271–280.50. Saltelli A. Making best use of model evaluations to compute sensitivity indices.Comp Phys Comm. 2002;145(2):280–297.51. Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M, Tarantola S. Variancebased sensitivity analysis of model output. Design and estimator for the totalsensitivity index. Comp Phys Comm. 2010;18(2):259–270.52. Aggarwal M, Hussaini MY, De La Fuente L, Navarrete F, Cogan NG. Aframework for model analysis across multiple experiment regimes: Investigatingeffects of zinc on
Xylella fastidiosa as a case study. Journal of TheoreticalBiology. 2018;457(14):88–100.53. Battista NA, Pearcy LB, Strickland WC. Modeling the prescription opioidepidemic. Bulletin of Mathematical Biology. 2019;81(7):2258–2289.doi:https://doi.org/10.1007/s11538-019-00605-0.54. Othmer HG, Dunbar SR, Alt W. Models of dispersal in biological systems.Journal of mathematical biology. 1988;26(3):263–298.
S1 Appendix Resource-Independence: The Telegrapher’s Equation.
One of our primary assertions is that, in order for the locust population to form acoherent traveling pulse, the two transition rates k sm , k ms must depend on the resourcedensity R . To demonstrate this we will estimate asymptotically the large-time behaviorfor 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 S t = − αS + βMM t = αS − βM − vM x x, t ∈ R (10)These equations are linear with constant coefficients and can be solved by a variety ofmeans. Physically, they correspond to the probability densities associated to anagent-based model of random switching between stationary and moving states. In thisinterpretation, α is the probability that an agent transitions from stationary to movingand β is the probability of transition from moving to stationary. As such, we canFebruary 26, 2020 30/37dentify this model as a variant of the Telegrapher’s Equation . Following previous work(see [54] and references therein), we use the method of moments to determine the largetime behavior.We take initial conditions corresponding to a starting at the origin, S ( x,
0) = S δ ( x ) , M ( x,
0) = M δ ( x ) (11)where S ( M ) is the probability is initially stationary (moving) and δ is the Dirac δ -function. We choose S + M = 1, reflecting the fact that this is a probability density.On average, an agent spends αα + β of its time moving. This leads us to conclude thatits average velocity is c ≡ v αα + β . This motivates a change of variables ξ = x − ct whichyields S t = − αS + βM + cS ξ M t = αS − βM + ( c − v ) M ξ x, t ∈ R (12)In this co-moving frame an agent now moves left with speed c (in the S state) and rightwith speed c − v (in the M state) but is never stationary.Solutions to this PDE correspond to probability distributions which can becharacterized by their moments. We define the n th moment M n ( t ) := (cid:90) ∞−∞ ξ n M ( ξ, t ) dξ (13) S n ( t ) := (cid:90) ∞−∞ ξ n S ( ξ, t ) dξ. (14)Multiplying (12) by ξ n and integrating yields the equations d S n dt = − α S n + β M n − nc S n − ( t ) d M n dt = α S n − β M n − n ( c − v ) M n − ( t ) t ∈ R + (15)A similar calculation yield new initial conditions S n (0) = S δ n, , M n (0) = M δ n, (16)where δ p,q is the Kronecker δ -function.When n = 0, we obtain the equations governing S and M the probability of abeing in state S (or M ) at time t , d S dt = − α S + β M , S (0) = S ,d M dt = α S − β M , M (0) = M . (17)The solution is S ( t ) = βα + β (cid:16) − e − ( α + β ) t (cid:17) + S e − ( α + β ) t (18) M ( t ) = αα + β (cid:16) − e − ( α + β ) t (cid:17) + M e − ( α + β ) t . (19)As t → ∞ , S ( t ) ∼ βα + β and M ( t ) ∼ αα + β as previously advertised.February 26, 2020 31/37he first moments S , M , when divided by S and M , are the centers of mass forthe corresponding the probability distributions S ( ξ, t ) , M ( ξ, t ). To find these, we solve d S dt = − α S + β M − c S ( t ) d M dt = α S − β M − n ( c − v ) M ( t ) t ∈ R + (20)and obtain a solution of the form S ( t ) = ( C + C t ) e − ( α + β ) t + βv ( α + β ) [ βM − α (1 + S )] (21) M ( t ) = ( C + C t ) e − ( α + β ) t + αv ( α + β ) [ − αS + β (1 + M )] (22)where the C i ’s are constants depending on v, α, β, S , M . We emphasize that each of S , M tend exponentially to a constant as t → ∞ . We conclude that, for large times,the centers of mass of the probability distributions are stationary in our co-movingframe. That is, they move with constant speed c = v αα + β as we hypothesized earlier.We repeat this procedure once more to find the second moments S , M , whichdescribe the variance of the probability distributions S ( ξ, t ) , M ( ξ, t ). The solution is ofthe form S ( t ) = ( D + D t + D t ) e − ( α + β ) t + D + (cid:18) αβ v ( α + β ) (cid:19) t (23) M ( t ) = ( D + D t + D t ) e − ( α + β ) t + D + (cid:18) α βv ( α + β ) (cid:19) t (24)where the D i ’s are again constants. Important here is that the variances of S ( ξ, t ) and M ( ξ, t ) are increasing linearly in t . This precludes the possibility of a coherent finitemass traveling pulse for large times.The Central Limit Theorem can be used to show that the distribution at large timesis normally distributed with linearly increasing variances. If we define the variances as( σ S ) = S ( t ) · S ( t ) − [ S ( t )] [ S ( t )] , ( σ M ) = M ( t ) · M ( t ) − [ M ( t )] [ M ( t )] then, as t → ∞ , ( σ S ) ∼ ( σ M ) ∼ ¯ σt, ¯ σ = 2 αβv ( α + β ) . and S ( ξ, t ) = β α + β ) √ π ¯ σt e − ξ σt + O ( t − / ) M ( ξ, t ) = α α + β ) √ π ¯ σt e − ξ σt + O ( t − / )which reinforces the random walk interpretation of this distribution.In summary, when k sm and k ms are independent of R , any initial densitydistribution will spread diffusively. Next, we will show that if dk sm dR ≤ dk ms dR ≥ S2 Appendix Traveling Wave Analysis.
Before analyzing the traveling wave solutions to the PDE (9), we simplify notationby rescaling variables into a nondimensional form. We apply the change of variables S = αλ ˜ S, M = αλ ˜ M , R = 1 γ ˜ R, x = vα ˜ x, t = 1 α ˜ t. (25)February 26, 2020 32/37ote that the dimensionless value of the initial boundary condition is ˜ R + = γR + andthe conserved quantity is now ˜ N = λα (cid:82) ∞−∞ S + M dx . This yields (after dropping the ∼ ’s) the dimensionless equations are R t = − SRS t = − k sm S + k ms MM t = k sm S − k ms M − M x (26)where now k sm = φ − ( φ − e − R and k ms = ψ − ( ψ − µ ) e − ωR (27)with four dimensionless parameters φ = ηα , ψ = θα , µ = βα , and ω = δγ . (28)Working from this non-dimensional PDE (26), we utilize a standard traveling-waveansatz. This amounts to assuming that solutions take the form S ( x, t ) = S ( x − ct ), andsimilarly for M, R . Here c is the speed of our moving reference frame for which we usethe spatio-temporal variable ξ = x − ct . The value of c is left to be determined in theforthcoming analysis. We obtain the ODEs R ξ = 1 c SRS ξ = 1 c ( k sm S − k ms M ) M ξ = 11 − c ( k sm S − k ms M ) . Subtracting the last two equations, we have cS ξ − (1 − c ) M ξ = 0. Integrating once, wehave the relation cS ( ξ ) − (1 − c ) M ( ξ ) = 0 (29)where we know the constant of integration on the right must be zero by considering thelong-time/far-distance limit ξ → ∞ . Using Eq (29), we rewrite the ODE above in termsof the total density ρ ( ξ ) = S ( ξ ) + M ( ξ ). We now have an equation amenable tophase-plane analysis R ξ = 1 − cc ρRρ ξ = (cid:18) k sm c − k ms − c (cid:19) ρ. (30)In this two-dimensional ODE, we prove the existence of heteroclinic connections thatcorrespond to traveling wave solutions of Eq (9). Theorem 1 (Existence of Traveling Waves) . For each c such that < φφ + ψ < c < µ < , there exists a one-parameter family of heteroclinic connectionsin the phase plane. We parameterize the family by R + . Each connection goes from from ( R − , to ( R + , and corresponds to a traveling wave solution to the PDE (9) thatmoves with speed c , and has uniquely determined total mass N , and leaves behindremaining resources R − . February 26, 2020 33/37 roof.
Consider the nullclines of system Eq (30). The R -nullclines are given by the lines R = 0 and ρ = 0. The ρ -nullclines are given by ρ = 0 and the vertical line that satisfies K ( R ) := k sm ( R ) c − k ms ( R )1 − c = φ − ( φ − e − R c − ψ − ( ψ − µ ) e − ωR − c = 0 . (31)The set of all equilibria is exactly the line ρ = 0; no other equilibria exist in the interiorof the first quadrant ρ, R > K (cid:48) ( R ) = 1 c dk sm dR − − c dk ms dR which is guaranteed to be less than zero for some c as long as dk sm dR ≤ dk ms dR ≥ ρ -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 + ). Notethat K (0) > K ( ∞ ) <
0, and K (cid:48) ( R ) <
0. By the continuity of P as R → ∞ , we canapply the Intermediate Value Theorem and guarantee a unique R ∗ ∈ (0 , R + ) for anylarge enough choice of R + .Fixing R + sufficiently large, we proceed with an invariant region argument, seeFigure 12. Define the rectangle A = { < R ≤ R ∗ , < ρ < ρ ∗ } for an arbitrary ρ ∗ to bedetermined below. Note that region A is invariant as ξ decreases. This is simply due tothe fact that R ξ > ρ ξ ≥ ρ -nullcline { R = R ∗ } must remain in region A as ξ decreases. By thePoincar´e-Bendixson Theorem, the trajectory must terminate on some point ( R − ,
0) as ξ → −∞ . A B
Fig 12.
The (
R, ρ ) phase plane with semi-invariant regions
A, B bounded by dottedlines including the ρ -nullcline (gold), sample arrows for the vector field (red), and aheteroclinic from ( R − ,
0) to ( R + ,
0) (blue).We define region B by choosing the upper boundary as a line segment parallel toand ε -above the stable eigenspace of the linearization of Eq (30) at ( R + ,
0) for someFebruary 26, 2020 34/37mall ε >
0. That is, B = { R ∗ ≤ R < R + , < ρ < (cid:96) ( R ) } where (cid:96) ( R ) = 1 R + c − c K ( R + )( R − R + ) + ε. (Note that we may now choose the upper bound of region A to be ρ ∗ = (cid:96) ( R ∗ ).)All that remains is to show that the stable manifold W s of ( R + ,
0) is containedwithin region B until it exits (as ξ decreases) through the vertical ρ -nullcline. Since W s is locally tangent to the stable eigenspace near ( R + , B as ξ → ∞ . Also W s cannot end, as ξ → −∞ , on any of the equilibriacomposing the lower boundary { ρ = 0 } of B . This is simply due to the fact that ρ ξ < B so we must have that ρ increases as ξ decreases. Since R ξ > R + ,
0) and ( R + , (cid:15) ), it is impossible for W s to exit B there as ξ decreases. Finally, W s cannot exit along (cid:96) because the slope of thevector field ρ ξ R ξ = 1 R c − c K ( R ) at any point ( R, ρ )is greater than the slope of (cid:96) . This follows from the fact that K (cid:48) ( R ) < R < R + on all of (cid:96) . In fact, we have shown that, as ξ decreases, no trajectory may exitregion 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 thatrelate
N, c, R + , and R − . Given any two of these variables and model parameter values,the following equations determine the other two variables: I I = cv − c (32) N λv = cv − c ln ( R + /R − ) (33)where I = (cid:90) R + R − k sm ( R ) R dR, and I = (cid:90) R + R − k ms ( R ) R dR. (34)We prove that equivalent equations hold for the nondimensionalized model Eq (26).
Theorem 2.
Given a traveling wave solution to Eq (26) . Then the quantities
N, c, R + , R − satisfy I I = c − c (35) N = c − c ln ( R + /R − ) (36) where I , I are given by Eq (34) with the nondimensional versions of k sm , k ms . Proof.
A traveling wave solution satisfies the ODE (30). Since R ξ > ρ as a function of R along any heteroclinic.Thus dρdξ = dρdR dRdξ so integrating along a heteroclinic, we have (cid:90) R + R − ρ ξ R ξ dR = (cid:90) R + R − dρdR dR = ρ ( R + ) − ρ ( R − ) = 0 . February 26, 2020 35/37e also have (cid:90) R + R − ρ ξ R ξ dR = (cid:90) R + R − c − c K ( R ) R dR = 11 − c I − c (1 − c ) I where K ( R ) is given in Eq (31). Therefore we have proved Eq (35).Dividing the equation for R ξ in Eq (30) by R and integrating the left hand side, wehave (cid:90) ∞−∞ R ξ R dξ = ln (cid:18) R + R − (cid:19) . Meanwhile, the right hand side gives us (cid:90) ∞−∞ − cc ρ dξ = 1 − cc N, proving Eq (36).In fact, these equalities can be used to prove monotonicity of the mass-speed relation. Theorem 3.
Fix R + . Then the speed c is a strictly increasing function of mass N anda strictly decreasing function of R − .Proof. Let s = c − c . Then Eqs (35)-(36) become s = I I (37) s = N ln (cid:16) R + R − (cid:17) . (38)First, we show dsdR − <
0. Taking the derivative of Eq (37), we get dsdR − = I dI dR − − I dI dR − I = I k ms ( R − ) − I k sm ( R − ) R − I . (39)We will show the numerator, call it I , is negative. Dividing I by k sm ( R − ) · k ms ( R − ),we get I k sm ( R − ) − I k ms ( R − ) = (cid:90) R + R − R (cid:18) k sm ( R ) k sm ( R − ) − k ms ( R ) k ms ( R − ) (cid:19) dR. (40)The integrand k sm ( R ) k s m ( R − ) − k m s ( R ) k m s ( R − ) is less than zero for R − < R < R + . To see this, notethat 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 negativesign) increases. Then the integrand is indeed negative for all R > R − .Second, we will show dsdN >
0. Differentiating Eqs (37) and (38) respectively, we get dSdN = ln (cid:16) R + R − (cid:17) + NR − dR − dN ln (cid:16) R + R − (cid:17) (41)and dSdN = dR − dN (cid:20) I k ms ( R − ) − I k sm ( R − ) I R − (cid:21) . (42)February 26, 2020 36/37etting these equal to each other, we obtain dR − dN = R − ln (cid:16) R + R − (cid:17) ln (cid:16) R + R − (cid:17) (cid:104) I k ms ( R − ) − I k sm ( R − ) I (cid:105) − N . (43)Substituting Eq (43) into either Eq (41) or Eq (42), we obtain dsdN = ln (cid:16) R + R − (cid:17) [ I k ms ( R − ) − I k sm ( R − )]ln (cid:16) R + R − (cid:17) [ I k ms ( R − ) − I k sm ( R − )] − N I = ln (cid:16) R + R − (cid:17) I ln (cid:16) R + R − (cid:17) I −
N I . (44)We have already shown that I <
0. Because R + > R − , the numerator is negative. Thedenominator is also negative, so have shown that dsdN >dsdN >