May 23, 2018, Wednesday
University of Colorado at Boulder Search A to Z Campus Map University of Colorado at Boulder CU 
Search Links

MBW:Locust Phase Change and Swarming

From MathBio

Jump to: navigation, search

Executive Summary

In Locust Dynamics: Behavioral Phase Change and Swarming[1], Topaz et. al discuss necessary conditions for an outbreak of marching hopper bands of locusts.
Locust swarm [2]
Outbreaks are triggered by a variety of factors, such as visual, olfactory, and tactile cues. Typically, locusts lead solitary lifestyles in arid regions. When overcrowding occurs, the locusts tend to change their behavior in the newly formed group to be gregarious, rather than solitarious. Topaz et. al developed a partial integrodifferential equation model in order to understand the dynamics of hopper band formation at the population level. The model is based on reaction-diffusion equations that relate the locusts’ behavioral phase changes and spatial movement. Stability analysis identifies conditions for an outbreak. Topaz et. al proceed to perform a model reduction in order to quantify the temporal dynamics, the subset of the population that will eventually group together, and the time until a swarm occurs. Finally, a numerical simulation reveals the structure of the swarm. They ultimately offer suggestions for locust control strategies:
  • Populations should be kept below a certain critical density threshold that may be calculated based on known biological parameters.
  • Prevention of the formation of gregarious bands is easier to achieve than controlling these swarms by means of extermination.
Artistic rendition of a swarm of desert locusts[3].

Context/Biological phenomenon under consideration

To see a locust swarm in action, check out this video.

Locusts are related to grasshoppers and crickets, but depending on their environment, they can change how they look and behave.
Adult locust [4]
Between outbreaks, in sparse surroundings, locusts live antisocial, solitary lifestyles. However, when resources are abundant, their numbers can easily rise to the point of overpopulation. This overcrowding induces morphological and behavioral changes in the locusts; they begin to tend towards social, aggregate behavior, their coloration becomes brighter, and their eyes, back legs, and wings enlarge to fly longer distances in search of food [2]. The authors succinctly describe a tactile cue that induces gregarious behavior in Schistocerca gregaria, the desert locust:

For the desert locust Schistocerca gregaria, the most potent stimulus is tactile: repetitive stroking of the femora of hind legs functions as a crowding indicator. Mechanosensory stimulation of leg nerves leads to serotonin cascades in the metathoracic ganglion, and initiates gregarious behavior. Gregarization can be induced by rubbing a locust’s hind leg for 5 s per minute during a period of 4 hr. Cessation of physical contact leads to solitarization after 4 hr.

The depleting resources in the area forces them to migrate en masse in search of new feeding grounds. This swarm of locusts can consist of millions of individuals, and it poses a huge threat to farmers. Swarms can travel a few hundred kilometers per day and span up to a thousand square kilometers [1]. Swarms are notorious for leaving farms without a single plant growing, and billions of dollars in crops are lost.


Previous locust models are based on Lagrangian simulations, where the position, velocity, and interactions of individual locusts are tracked. Also, the focal point of many current models is interactions with clusters of resources and varying environmental conditions. Furthermore, some models use anisotropic interactions (different responses to anterior or posterior neighbors) or consider Newtonian dynamics[1]. Topaz et. al present a model that is based on population density (Eulerian), rather than individuals. This allows them to use techniques of partial differential equations and their extensions (integro-PDEs). They consider the attractive-repulsive social interactions between locusts in their model instead of interactions between resource clusters and changing environmental conditions. Lastly, the authors rely on isotropic interactions (as opposed to anisotropic ones) in their model.


  • The model uses a density field of locusts rather than tracking individuals.
  • Individuals can sense density nearby in terms of attractive-repulsive forces and social potentials.
  • Social interaction is radially symmetric.
  • Nonlocal interaction models capture spatially distributed interactions.
  • The model is spatially periodic.

Model construction

Classic swarm modeling involves a conserved population density field \rho travelling at a velocity {\mathbf  {v}} that arises from social interactions:

\rho _{t}+\nabla \cdot (\rho {\mathbf  {v}})=0

The contribution of \rho ({\mathbf  {x'}},t) of a small group of individuals at location {\mathbf  {x'}} to the force on the individual at position {\mathbf  {x}} is given by:

{\mathbf  {F(x-x')}}\rho ({\mathbf  {x'}},t)=-\nabla Q{\mathbf  {(x-x')}}\rho ({\mathbf  {x'}},t)


{\mathbf  {v(x,}}t)=-\int _{\Omega }\nabla Q{\mathbf  {(x-x')}}\rho ({\mathbf  {x'}},t)d{\mathbf  {x'}}

The reader should note that the term Q{\mathbf  {(x-x')}} is radially symmetric. Let us use some notation to represent this integral:

{\mathbf  {v}}=-\nabla Q*\rho

This is a nonlocal interaction model, and its purpose is to capture spatially distributed interactions, as opposed to pure PDEs (which describe interactions over infinitesimal ranges).

This basic model is further adapted for biphasic insects below (the Full model).


  • \rho (x,t) is the density field of locusts
  • F attractive-repulsive forces
  • Q social potentials; these can be solitarious or gregarious
  • \Omega represents the spatial domain
  • {\mathbf  {x}}=(x,y) is the spatial coordinate

Full model


  • Locust sensing Q_{{s,g}} is directionally isotropic (uniform in all directions).
  • Repulsion dominates at small length scales.
  • Attraction dominates at larger length scales.
  • The environment is spatially homogeneous.
  • There is no reproductive or death term in the equations because the time scale for the formation of a swarm is relatively short compared to the life expectancy of a locust and its reproductive time scale.

Model construction

There are separate density fields for the solitarious s({\mathbf  {x}},t) and gregarious phases g({\mathbf  {x}},t), and these phases sum to the total density:

\rho =s({\mathbf  {x}},t)+g({\mathbf  {x}},t)

The reaction-diffusion equations of the model are:

{\dot  {s}}+\nabla \cdot ({\mathbf  {v_{s}}}s)=-f_{2}(\rho )s+f_{1}(\rho )g

{\dot  {g}}+\nabla \cdot ({\mathbf  {v_{g}}}g)=f_{2}(\rho )s-f_{1}(\rho )g

Where the velocities are:

{\mathbf  {v_{s}}}=-\nabla (Q_{s}*\rho )

{\mathbf  {v_{g}}}=-\nabla (Q_{g}*\rho )

We now require equations for the solitarious/gregarious social interactions Q_{{s,g}} and the two density-dependent conversion rate functions, f_{{1,2}}(\rho ), which represent the rate of gregarious \rightarrow solitarious and the rate of solitarious \rightarrow gregarious, respectively:

Q_{s}({\mathbf  {x-x'}})=R_{s}e^{{-|{\mathbf  {x-x'}}|/r_{s}}}

Failed to parse (lexing error): Q_g(\mathbf{ x-x'}) = R_g e^{-|\mathbf{x-x'}|/r_g} –- A_g e^{-|\mathbf{x-x'}|/a_g}

f_{1}(\rho )={\frac  {\delta _{1}}{1+(\rho /k_{1})^{2}}}

f_{2}(\rho )={\frac  {\delta _{2}(\rho /k_{2})^{2}}{1+(\rho /k_{2})^{2}}}

Conversion rates from one phase to the other

By biological observation, we know that the time to gregarization is relatively fast compared to the time to solitarization. Thus, the equations for f_{{1,2}} are justified. It should be noted that f_{1} decreases with \rho , while f_{2} increases with \rho to the saturation point, \delta _{2}.

Critical distance

As indicated above, it should be noted that Q_{s} is a decaying exponential function for all choices of R_{s} and s_{s}, so it is purely repulsive. However, Q_{g}, being the difference of two exponentials, suggests that there is an equilibrium, where there is no net contribution to the velocity. The equilibrium point d is found by solving -\nabla Q_{g}=0.

d={\frac  {a_{g}r_{g}}{a_{g}-r_{g}}}\ln({\frac  {R_{g}a_{g}}{A_{g}r_{g}}})

d is also called the critical distance. d refers to the spatial locations {\mathbf  {x',x}}, so it does not describe population-level features. It is important to note that d will be negative for certain values of R_{g},A_{g},r_{g},a_{g}, which is physically impossible. Furthermore, even for positive values of d, individuals may disperse, aggregate, or clump unexpectedly. We must choose conditions such that the tendency of locusts to aggregate is modeled accurately. Therefore, we impose these restrictions on R_{g},A_{g},r_{g},a_{g}:

For repulsion to dominate at small length scales,

Failed to parse (lexing error): R_ga_g –- A_gr_g > 0

For attraction to dominate at larger length scales,

Failed to parse (lexing error): A_ga_g^2 –- R_gr_g^2 > 0

These conditions are based on existing literature that shows cohesiveness requires parameters in Q_{g} must be such that locusts are led to clumping [1].


The numerical values of these parameters are based on estimates from biological literature [1].

  • \rho is total density
  • s({\mathbf  {x}},t) is the density field for solitarious individuals
  • g({\mathbf  {x}},t) is the density field for gregarious individuals
  • f_{1} rate of phase change from gregarious \rightarrow solitarious
  • f_{2} rate of phase change from solitarious \rightarrow gregarious
  • \delta _{{1,2}} maximal phase transition rates
    • For situations where rates of gregarization and solitarization occur on the same time scale: \delta _{{1,2}}=\delta =0.25hr^{{-1}}
    • For situations with large differences in rates of gregarization and solitarization: \delta _{1}=0.025hr^{{-1}},\delta _{2}=0.25hr^{{-1}}
  • k_{{1,2}} characteristic locust densities where f_{{1,2}} are half of their maximal values
    • For situations where the critical density for the onset of collective motion is 50-80locusts/m^{2}: k_{{1,2}}=k=65locusts/m^{2}
    • For situations where the gregarious \rightarrow solitarious transition occurs at a higher density threshold than the solitarious \rightarrow gregarious transition: k_{1}=20locusts/m^{2},k_{2}=65locusts/m^{2}
  • Q_{s} solitarious social interaction, purely repulsive
  • Q_{g} gregarious social interaction
  • R_{s},R_{g},A_{g} interaction amplitudes that determine the strengths of attraction and repulsion
    • R_{s}=11.87m^{3}/(hr\cdot locust)
    • R_{g}=5.13m^{3}/(hr\cdot locust)
    • A_{s}=13.33m^{3}/(hr\cdot locust)
  • r_{s},r_{g},a_{g} interaction length scales that represent typical distances over which one locust can sense and respond to another
    • r_{g}=0.04m (repulsion length scale for gregarious individuals)
    • r_{s}=0.14m (repulsion length scale for solitarious individuals)
    • a_{g}=0.14m (attraction length scale)


Topaz et. al fist find the homogeneous steady state solution to the reaction-diffusion equations and analyze the solutions using linear stability analysis, which is equivalent to finding the eigenvalues of the linearized system. The stability analysis reveals the condition for the onset of a swarm. The authors then proceed to perform numerical simulations in one spatial dimension to identify the structure of the swarm. They find the model exhibits hysteresis on the population level (macroscopic properties of the reaction-diffusion equations, which are outputs of the model).

Homogeneous steady states (HSS)

The steady state solutions to the reaction-diffusion equations for low density \rho _{0} are:

s_{0}\approx \rho _{0}-{\frac  {\delta _{2}}{\delta _{1}k_{2}^{2}}}\rho _{0}^{3}

g_{0}\approx {\frac  {\delta _{2}}{\delta _{1}k_{2}^{2}}}\rho _{0}^{3}

For large density \rho _{0}:

s_{0}\approx {\frac  {\delta _{1}k_{1}^{2}}{\delta _{2}\rho _{0}}}

g_{0}\approx \rho _{0}-{\frac  {\delta _{1}k_{1}^{2}}{\delta _{2}\rho _{0}}}

Using the default parameters k=65,\delta =0.25, the following plot A is a graph of the steady state solutions of s_{0},g_{0} as total density \rho _{0} is varied, in blue and green, respectively. Plot B is a zoomed in view of the region near \rho _{*}, where g_{0} becomes larger than s_{0}:


The middle curve in each set of curves is the solution for the default choice of parameters. The other curves show parameter sensitivity. In the gray and red regions, the HSS is linearly unstable to small perturbations. At the critical density \rho _{*}, s_{0} reaches a maximum, and g_{0} continues to grow. As a result of the authors’ choice of the default parameters, s_{0}^{{max}}=k/2 occurs at \rho _{*}=k, but this is not always the case, as is seen in the top and bottom curves (the 25th and 75th percentile values).

The general case where \delta _{1}\neq \delta _{2},k_{1}\neq k_{2} is expressed below, where the sum of the fractions of solitarious and gregarious locusts equals one:

\phi _{s}+\phi _{g}=1

The gregarious state necessarily becomes the dominant one as density increases:

\phi _{g}={\frac  {1}{1+\gamma K^{2}{\frac  {1+\psi ^{2}}{\psi ^{2}(\psi ^{2}+K^{2})}}}}

  • \gamma ={\frac  {\delta _{1}}{\delta _{2}}} ratio of maximal solitarization rate to maximal gregarization rate
  • K={\frac  {k_{1}}{k_{2}}} ratio of the characteristic solitarization and gregarization densities
  • \psi ={\frac  {\rho _{0}}{k_{2}}} rescaled spatially homogeneous density

Linear stability analysis

Linear stability results depend on Fourier transforms {\hat  {Q}}_{{s,g}}(q) of the interaction potentials Q_{{s,g}}. The eigenvalue \lambda determines the stability of the locust model:

{\hat  {Q}}_{s}(q)={\frac  {2\pi R_{s}r_{s}^{2}}{(1+r_{s}^{2}q^{2})^{{3/2}}}}

{\hat  {Q}}_{g}(q)={\frac  {2\pi R_{g}r_{g}^{2}}{(1+r_{g}^{2}q^{2})^{{3/2}}}}-{\frac  {2\pi A_{g}a_{g}^{2}}{(1+a_{g}^{2}q^{2})^{{3/2}}}}

\lambda _{1}(q)=-q^{2}[s_{0}{\hat  {Q}}_{s}(q)+g_{0}{\hat  {Q}}_{g}(q)]

  • q=|{\mathbf  {q}}| perturbation wave number

Based on the equations above, Topaz et. al found the instability condition in terms of \phi _{g}:

Failed to parse (lexing error): \phi_g > \phi_g^*=\frac{R_sr_s^2}{R_sr_s^2 – R_gr_g^2 +A_ga_g^2 }

Therefore, if a sufficiently large proportion of the population is gregarious, the solution is unstable. Substitution of \phi _{g}^{*} into the general solution of \phi _{g} above will yield a solution in terms of \rho _{o}. That algebra is tedious, so the authors refer to a contour plot instead. Due to the fact that both \phi _{g},\rho _{0} are monotonically related, the HSS is unstable for sufficiently dense populations (large values of \rho _{0}). The contour plot below illustrates the stability of the HSS:


\gamma ,K vary along the horizontal and vertical axes. The contours are critical values of \psi ^{*}, the rescaled density. For \psi >\psi ^{*}, the solutions are unstable. The arrow along the horizontal axis shows where \gamma moves as the rate of gregarization increases. The arrow along the vertical axis shows where K travels as the density threshold for gregarization decreases. It should be noted that one should use \rho ^{*} for an accurate biological interpretation, not \psi ^{*} (\rho ^{*}=\psi ^{*}k_{2}).

The default parameters yield the following results:

  • The homogeneous solution is unstable for \rho >\rho ^{*}=62.3locusts/m^{2} (the value in the gray region of Plot B).
  • \phi _{g}^{*}=0.479
  • K=\gamma =1
  • \psi =0.959

The authors found extremely disparate results with different parameters, K=1/3,\gamma =1/10:

  • \rho _{0}^{*}=15.9locusts/m^{2}

Using the default parameters (blue and dashed green curves of plots A, B), they found the wave number q_{{max}} corresponding to the most rapidly growing perturbation. Their findings are shown in the plot below:


Again, the middle q_{{max}} curve of this plot matches the middle curve in plots A,B. For the middle curve, saturation occurs at q_{{max}}=8.89m^{{-1}}. This value indicates fast-growing perturbations occur on the length scale of a few locust bodies because the length scale 2\pi /q_{{max}}\approx 0.71m.

Unfortunately, this linear stability analysis only applies to the case where we observe small perturbations of uniform steady states. It does not predict long-term or large-amplitude dynamics, and neither does it accurately describe large perturbations. Furthermore, using this model to analyze small perturbations in a nonuniform steady state would also be erroneous.

Numerical simulation

Topaz et. al construct a one dimensional numerical simulation to examine the structure of the swarm. The parameters of the simulation are as follows:

  • L=3m domain length (domain is periodic)
  • M=50locusts total population mass
  • R_{s}=6.83m^{2}/(hr\cdot locust)
  • R_{g}=6.04m^{2}/(hr\cdot locust)
  • A_{g}=12.9m^{2}/(hr\cdot locust)
  • r_{s},r_{g},a_{g} are the same as before
  • For the default parameters:
    • k_{{1,2}}=k=8locusts/m
    • \delta _{{1,2}}=\delta =0.25hr^{{-1}}
  • For the alternative parameters:
    • k_{1}=4.5locusts/m,k_{2}=8locusts/m
    • \delta _{1}=0.025hr^{{-1}},\delta _{2}=0.25hr^{{-1}}

The results are shown in the figures below, where the first one is with the default set of parameters and the second one is with the alternative set of parameters:



This snapshots of the number of solitarious and gregarious locusts at different times illustrates the transition from a population of initially solituarious locusts (t=0hr) to a gregarious swarm (t=10,80hr).


The numerical results suggest the model has population-level hysteresis for the alternative parameters (macroscopic properties of the reaction-diffusion equations, which are outputs of the model):


Unstable behavior is represented by dashed lines and stable behavior is represented by the solid lines. The red curve is the HSS solution, where there is no clustering. The green curve represents the clustered state from the numerical simulation. The asterisk on the horizontal axis represents the point of linear instability, where stability jumps from the red curve to the green curve.

As \rho _{0} increases, the stable solution is along the red curve. As \rho _{0} passes through the asterisk, stability changes to the green curve. As \rho _{0} is decreased from this point, stability remains on the green curve (even for small values of \rho _{0}).

Hysteresis suggests that once a locust swarm has formed, it cannot be eliminated by reducing the density of the swarm.


The results of this investigation can be summarized as follows:

  • Locusts exist in a spatially uniform steady state where the density of locusts is less than a critical density. Beyond this point, the stable distribution is dense gregarious clusters (swarms).
  • Stability analysis shows that this critical density depends on ratios of biological parameters.
  • Numerical simulations show how gregarization progresses within a matter of hours.
  • Eventually, the entire population becomes gregarious.
  • Population-level hysteresis shows that past the critical density, where the distribution is stable in clusters, even great reductions of population density do not shift locust behavior from gregarious to solitarious.

Topaz et. al suggest the following locust management strategies:

  • Locust populations should be kept below the critical density in order to avoid an outbreak of swarming.
  • Exterminating a gregarious subset of individuals is more difficult than preventing group formation; prevention is easier than control in this case.

External links


  1. 1.0 1.1 1.2 1.3 1.4 Topaz, C., D’Orsongna, M., Edelstein-Keshet, L., Bernoff, A. (2012). Locust Dynamics: Behavioral Phase Change and Swarming. PLOS Computational Biology, 8.
  2. 2.0 2.1