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

MBW:Two Models of Locust Phase Change

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. 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 solitary. This article explores a modification of Topaz et. al's original reaction-diffusion model of the behavior phase change dynamics. Stability analysis of the original model and the modified model identifies conditions for an outbreak. Specifically, this investigation explores the case where functions that describe the rate of change of phase conversion are modified to occur at faster rates (lower population densities). The main findings are:

  • Populations should be kept below a certain critical density threshold that may be calculated based on known biological parameters.
  • Choosing various functions that model phase change from solitary to gregarious behavior is not as influential in stability analysis as choosing various social interaction potential functions.
A desert locust[2].

Context/Biological phenomenon under consideration

Locusts are related to grasshoppers and crickets, but depending on their environment, they can change how they look and behave. Between outbreaks, locusts live antisocial, solitary lifestyles in sparse surroundings with low population densities. However, when there is an abundance of resources, populations rise to the point of overpopulation. Increasing population densities induce 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. The depleting resources in the area forces them to migrate en masse in search of new feeding and breeding 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. Interestingly, locusts at the rear of the swarm are known to exhibit cannibalistic behavior as locusts at the front of the swarm devour all the vegetation in the area. Without anything left to eat, the locusts at the rear consume the locusts directly in front of them as the swarm travels forwards. An interesting article from ScienceDaily about cannibalistic tendencies in locusts is worth reading.

The following video shows how a swarm of locusts travel.

Alt text
Locust life cycle [3]


For more information about the history of previous locust swarming models, consult the wiki page MBW:Locust Phase Change and Swarming. The model Topaz et. al explore is based on population density (rather than individuals), social interaction potentials between locusts, and expressions that describe the rate a solitary locust becomes gregarious and vice versa. These qualities allow them to use techniques of PDEs and their extensions (integro-PDEs). Below, the classic swarm model is described. The findings of this investigation are based on a modified version of this 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 in terms of the environment the locust swarm.

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).


  • \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}}}

Alt text
Plot of f_{{1,2}}(\rho )

The modified version of the model preserves the same expressions to describe velocity and the social interaction potential. Where this model differs from the original model in Topaz et. al is in the conversion rate functions, f_{{1,2}}. The modified functions are:

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

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

Alt text
Plot of {\tilde  {f}}_{{1,2}}(\rho )
Alt text
Plot of {\tilde  {f}}_{{1,2}}(\rho ) (zoomed in)

The purpose of introducing exponential functions {\tilde  {f}}_{{1,2}} is to determine whether the rate locusts change behavior affects the stability of the overall model. These exponentially growing and decaying functions model locust behavior change at a faster rate than the ones from the original paper. This investigation determines the extent to which behavior conversion rates affect the system's overall stability. It should be noted that both f_{{1,2}} and {\tilde  {f}}_{{1,2}} are monotonically decreasing and increasing, respectively.

Critical distance

Finding the critical distance leads to certain requirements in the model with respect to the interaction amplitudes and interaction length scales, R_{{g,s}},a_{g},A_{g},r_{{g,s}}. These calculations are based on the social interaction potentials, Q_{{s,g}}. As the modified model preserves the original social interaction potentials, the original analysis from Topaz et. al regarding critical distance remain valid.

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


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

  • \rho is total density
  • s,g is the density field for solitary/gregarious individuals
  • f_{1},{\tilde  {f}}_{1} rate of phase change from gregarious \rightarrow solitarious
  • f_{2},{\tilde  {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,g}} solitary/gregarious social interaction
  • R_{s},R_{g},A_{g} interaction amplitudes that determine the strengths of attraction and repulsion
  • r_{s},r_{g},a_{g} interaction length scales that represent typical distances over which one locust can sense and respond to another


The analysis of these two models has two parts:

  • Homogeneous steady states (HSS)
  • Linear stability

Homogeneous steady states

The steady states of the model can be obtained by setting the derivatives in time and space to zero and using the expression for \rho .

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

0=-f_{2}(\rho )s+f_{1}(\rho )g

0=f_{2}(\rho )s-f_{1}(\rho )g

Thus, the expressions for the equilibrium points s_{0},g_{0} are:

s_{0}=g_{0}{\frac  {f_{1}}{f_{2}}}

g_{0}={\frac  {\rho }{{\frac  {f_{1}}{f_{2}}}+1}}

For the original model, the full expression is:

s_{0}={\frac  {\rho _{0}\delta _{1}k_{1}^{2}(k_{2}^{2}+\rho _{0}^{2})}{\delta _{1}k_{1}^{2}(k_{2}^{2}+\rho _{0}^{2})+\delta _{2}rho_{0}^{2}(k_{1}^{2}+\rho _{0}^{2})}}

g_{0}={\frac  {\rho _{0}^{3}\delta _{2}(k_{1}^{2}+\rho _{0}^{2})}{\delta _{1}k_{1}^{2}(k_{2}^{2}+\rho _{0}^{2})+\delta _{2}rho_{0}^{2}(k_{1}^{2}+\rho _{0}^{2})}}

For the case where k_{1}=k_{2}, this is a plot of the solutions:

Alt text
HSS solutions of the original model, where some parameters are equal. The intersection of the two curves occurs at \rho _{0}\approx 65

For the case where k_{1}\neq k_{2}, this is a plot of the solutions:

Alt text
HSS solutions of the original model, where the parameters are distinct. The intersection of the two curves occurs at \rho _{0}\approx 8

The corresponding equations for the modified model:

{\tilde  {s}}_{0}={\frac  {\rho _{0}\delta _{1}(1+e^{{\rho _{0}/k_{2}}})}{\delta _{1}(1+e^{{\rho _{0}/k_{2}}})+\delta _{2}(e^{{{\frac  {\rho _{0}}{k_{1}}}+{\frac  {2\rho _{0}}{k_{2}}}}}+e^{{2\rho _{0}/k_{2}}})}}

{\tilde  {g}}_{0}={\frac  {\rho _{0}\delta _{2}e^{{2\rho _{0}/k_{2}}}(1+e^{{\rho _{0}/k_{1}}})}{\delta _{1}(1+e^{{\rho _{0}/k_{2}}})+\delta _{2}(e^{{{\frac  {\rho _{0}}{k_{1}}}+{\frac  {2\rho _{0}}{k_{2}}}}}+e^{{2\rho _{0}/k_{2}}})}}

For the case where k_{1}=k_{2}, this is a plot of the solutions:

Alt text
HSS solutions of the modified model, where some parameters are equal.

The following is a zoomed in plot near the origin:

Alt text
HSS solutions of the modified model, where some parameters are equal (zoomed in).

For the case where k_{1}\neq k_{2}, this is a plot of the solutions:

Alt text
HSS solutions of the modified model, where the parameters are distinct.

The following is a zoomed in plot near the origin:

Alt text
HSS solutions of the modified model, where the parameters are distinct (zoomed in).

Linear stability

The process of performing a linear stability analysis on the HSS solutions is outlined below:

  1. Consider small perturbations
  2. Linearize
  3. Fourier expansion
  4. Construct matrix of each Fourier mode amplitude
  5. Eigenvalues determine instability conditions

Consider small perturbations

The perturbations about the HSS solutions s_{0},g_{0} are denoted as s_{1},g_{1}:

s({\mathbf  {x}},t)=s_{0}+s_{1}({\mathbf  {x}},t)

g({\mathbf  {x}},t)=g_{0}+g_{1}({\mathbf  {x}},t)

Accordingly, our expression for \rho _{0} changes slightly.

\rho _{0}({\mathbf  {x}},t)=s_{0}+s_{1}({\mathbf  {x}},t)+g_{0}+g_{1}({\mathbf  {x}},t)

For convenience, refer to this table that shows the analysis of the original model and the modified model side by side:

Original model Modified model
s({\mathbf  {x}},t)=s_{0}+s_{1}({\mathbf  {x}},t) {\tilde  {s}}({\mathbf  {x}},t)={\tilde  {s}}_{0}+{\tilde  {s}}_{1}({\mathbf  {x}},t)
g({\mathbf  {x}},t)=g_{0}+g_{1}({\mathbf  {x}},t) {\tilde  {g}}({\mathbf  {x}},t)={\tilde  {g}}_{0}+{\tilde  {g}}_{1}({\mathbf  {x}},t)
\rho _{0}({\mathbf  {x}},t)=s_{0}+s_{1}({\mathbf  {x}},t)+g_{0}+g_{1}({\mathbf  {x}},t) {\tilde  {\rho }}({\mathbf  {x}},t)={\tilde  {s}}_{0}+{\tilde  {s}}_{1}({\mathbf  {x}},t)+{\tilde  {g}}_{0}+{\tilde  {g}}_{1}({\mathbf  {x}},t)


Substitution and further simplification of the perturbed equations into the reaction-diffusion equations leads to a linearized system. For reference, the reaction diffusion equations 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

Original model Modified model
{\dot  {s}}=s_{0}Q_{s}*\nabla ^{2}(s_{1}+g_{1})-As_{1}+Bg_{1} {\dot  {{\tilde  {s}}}}={\tilde  {s}}_{0}Q_{s}*\nabla ^{2}({\tilde  {s}}_{1}+{\tilde  {g}}_{1})-{\tilde  {A}}{\tilde  {s}}_{1}+{\tilde  {B}}{\tilde  {g}}_{1}
{\dot  {g}}=g_{0}Q_{g}*\nabla ^{2}(s_{1}+g_{1})+As_{1}-Bg_{1} {\dot  {{\tilde  {g}}}}={\tilde  {g}}_{0}Q_{g}*\nabla ^{2}({\tilde  {s}}_{1}+{\tilde  {g}}_{1})+{\tilde  {A}}{\tilde  {s}}_{1}-{\tilde  {B}}{\tilde  {g}}_{1}
A=f_{2}+f_{2}'s_{0}-f_{1}'g_{0} {\tilde  {A}}={\tilde  {f}}_{2}+{\tilde  {f}}_{2}'{\tilde  {s}}_{0}-{\tilde  {f}}_{1}'{\tilde  {g}}_{0}
B=f_{1}+f_{1}'g_{0}-f_{2}'s_{0} {\tilde  {B}}={\tilde  {f}}_{1}+{\tilde  {f}}_{1}'{\tilde  {g}}_{0}-{\tilde  {f}}_{2}'{\tilde  {s}}_{0}

Note that Q_{{s,g}} is present in both the original and the modified model. f_{{1,2}},{\tilde  {f}}_{{1,2}} are monotonically decreasing and increasing, respectively. It follows that A,{\tilde  {A}},B,{\tilde  {B}}>0 for all \rho _{0},{\tilde  {\rho }}_{0}>0.

Fourier expansion

A Fourier expansion is required to further analyze the linear system. We expand the perturbations as an infinite Fourier series. The domain is infinitely large. {\mathbf  {q}} is the wave number.

Original model Modified model
s_{1}({\mathbf  {x}},t)=\sum \limits _{{{\mathbf  {q}}}}S_{{\mathbf  {q}}}(t)e^{{i{\mathbf  {qx}}}} {\tilde  {s}}_{1}({\mathbf  {x}},t)=\sum \limits _{{{\mathbf  {q}}}}{\tilde  {S}}_{{\mathbf  {q}}}(t)e^{{i{\mathbf  {qx}}}}
g_{1}({\mathbf  {x}},t)=\sum \limits _{{{\mathbf  {q}}}}G_{{\mathbf  {q}}}(t)e^{{i{\mathbf  {qx}}}} {\tilde  {g}}_{1}({\mathbf  {x}},t)=\sum \limits _{{{\mathbf  {q}}}}{\tilde  {S}}_{{\mathbf  {q}}}(t)e^{{i{\mathbf  {qx}}}}

Construct matrix of each Fourier mode amplitude

Using the Fourier expansions in the linearized system reveals a ODEs for each Fourier mode amplitude. The Fourier transforms of the social interaction potentials, Q_{{s,g}} are defined as:

{\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}}}}

The social interaction potentials and their respective Fourier transforms are the same for both the original model and the modified model.

The matrix form of our linear system:

Original model Modified model
{\frac  {d}{dt}}\left({\begin{matrix}S_{{\mathbf  {q}}}\\G_{{\mathbf  {q}}}\end{matrix}}\right)={\mathbf  {L}}(q)\left({\begin{matrix}S_{{\mathbf  {q}}}\\G_{{\mathbf  {q}}}\end{matrix}}\right) {\frac  {d}{dt}}\left({\begin{matrix}{\tilde  {S}}_{{\mathbf  {q}}}\\{\tilde  {G}}_{{\mathbf  {q}}}\end{matrix}}\right)={\tilde  {{\mathbf  {L}}}}(q)\left({\begin{matrix}{\tilde  {S}}_{{\mathbf  {q}}}\\{\tilde  {G}}_{{\mathbf  {q}}}\end{matrix}}\right)
{\mathbf  {L}}(q)=\left({\begin{matrix}-s_{0}q^{2}{\hat  {Q}}_{s}(q)-A&-s_{0}q^{2}{\hat  {Q}}_{s}(q)+B\\-g_{0}q^{2}{\hat  {Q}}_{g}(q)+A&-g_{0}q^{2}{\hat  {Q}}_{g}(q)-B\end{matrix}}\right) {\tilde  {{\mathbf  {L}}}}(q)=\left({\begin{matrix}-{\tilde  {s}}_{0}q^{2}{\hat  {Q}}_{s}(q)-{\tilde  {A}}&-{\tilde  {s}}_{0}q^{2}{\hat  {Q}}_{s}(q)+{\tilde  {B}}\\-{\tilde  {g}}_{0}q^{2}{\hat  {Q}}_{g}(q)+{\tilde  {A}}&-{\tilde  {g}}_{0}q^{2}{\hat  {Q}}_{g}(q)-{\tilde  {B}}\end{matrix}}\right)

Eigenvalues determine instability conditions

The eigenvalues of the matrix reveal stability conditions for the system. Positive eigenvalues indicate growing perturbations, which lead to instability. Negative eigenvalues correspond with stable systems.

Original model Modified model
\lambda _{1}(q)=-q^{2}[s_{0}{\hat  {Q}}_{s}(q)+g_{0}{\hat  {Q}}_{g}(q)] {\tilde  {\lambda }}_{1}(q)=-q^{2}[{\tilde  {s}}_{0}{\hat  {Q}}_{s}(q)+{\tilde  {g}}_{0}{\hat  {Q}}_{g}(q)]
\lambda _{2}(q)=-(A+B) {\tilde  {\lambda }}_{2}(q)=-({\tilde  {A}}+{\tilde  {B}})

Since \lambda _{2},{\tilde  {\lambda }}_{2} are always negative, we look to \lambda _{1},{\tilde  {\lambda }}_{1} for the instability conditions. In other words, we are looking for the point where \lambda _{1},{\tilde  {\lambda }}_{1}>0.

At this point, let us rewrite s_{0},g_{0},{\tilde  {s}}_{0},{\tilde  {g}}_{0} as mass fractions that sum to 1 (for convenience).

Original model Modified model
\phi _{g}+\phi _{s}=1 {\tilde  {\phi }}_{g}+{\tilde  {\phi }}_{s}=1
\phi _{g}={\frac  {g_{0}}{g_{0}+s_{0}}} {\tilde  {\phi }}_{g}={\frac  {{\tilde  {g}}_{0}}{{\tilde  {g}}_{0}+{\tilde  {s}}_{0}}}

\lambda _{1},{\tilde  {\lambda }}_{1} now become:

Original model Modified model
\lambda _{1}(q)=-q^{2}[(1-\phi _{g}){\hat  {Q}}_{s}(q)+\phi _{g}{\hat  {Q}}_{g}(q)] {\tilde  {\lambda }}_{1}(q)=-q^{2}[(1-{\tilde  {\phi }}_{g})_{0}{\hat  {Q}}_{s}(q)+{\tilde  {\phi }}_{g}{\hat  {Q}}_{g}(q)]

The point where \lambda _{1},{\tilde  {\lambda }}_{1}>0 change sign occurs where the first term in the square brackets becomes smaller than the second term. Expansion and rearranging the eigenvalues by factoring out the locust attraction term from {\hat  {Q}}_{g} of the expression

Original model Modified model
\phi _{g}{\frac  {2\pi A_{g}a_{g}^{2}}{(1+a_{g}^{2}q^{2})^{{3/2}}}} {\tilde  {\phi }}_{g}{\frac  {2\pi A_{g}a_{g}^{2}}{(1+a_{g}^{2}q^{2})^{{3/2}}}}

leaves us with:

\lambda _{1}(q)=-\rho _{0}q^{2}\phi _{g}{\frac  {2\pi A_{g}a_{g}^{2}}{(1+a_{g}^{2}q^{2})^{{3/2}}}}[{\frac  {1-\phi _{g}}{\phi _{g}}}{\frac  {R_{s}r_{s}^{2}(1+a_{g}^{2}q^{2})^{{3/2}}}{A_{g}a_{g}^{2}(1+r_{s}^{2}q^{2})^{{3/2}}}}+{\frac  {R_{g}r_{g}^{2}(1+a_{g}^{2}q^{2})^{{3/2}}}{A_{g}a_{g}^{2}(1+r_{g}^{2}q^{2})^{{3/2}}}}-1]

{\tilde  {\lambda }}_{1}(q)=-\rho _{0}q^{2}{\tilde  {\phi }}_{g}{\frac  {2\pi A_{g}a_{g}^{2}}{(1+a_{g}^{2}q^{2})^{{3/2}}}}[{\frac  {1-{\tilde  {\phi }}_{g}}{{\tilde  {\phi }}_{g}}}{\frac  {R_{s}r_{s}^{2}(1+a_{g}^{2}q^{2})^{{3/2}}}{A_{g}a_{g}^{2}(1+r_{s}^{2}q^{2})^{{3/2}}}}+{\frac  {R_{g}r_{g}^{2}(1+a_{g}^{2}q^{2})^{{3/2}}}{A_{g}a_{g}^{2}(1+r_{g}^{2}q^{2})^{{3/2}}}}-1]

The expressions for \lambda _{1},{\tilde  {\lambda }}_{1} are essentially identical at this point, excepting the \phi _{g},{\tilde  {\phi }}_{g} terms. Furthermore, the two functions increase monotonically with \rho _{0}, so whatever the implications in terms of density are for one function, the same implications for the other follow. Thus, from this point onwards we will proceed in terms of \phi _{g}. For the first term in the square brackets to be small, the factor of {\frac  {1-\phi _{g}}{\phi _{g}}} must be small. This implies that \phi _{g} should be large, which correlates with increasing density.

If the bracketed quantity should change sign, it will happen where the first two terms are minimized, since these are positive terms. It is biologically reasonable to assume that a_{g}>=r_{s} and a_{g}>r_{g}. Therefore, the first two terms of the bracketed quantity are either constant or monotonically increasing. Therefore, the minimum occurs where q=0. Using q=0 in the expression for \lambda _{1} leaves us with:

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

Rearranging this expression in terms of density would necessarily lead to different expressions for \rho _{0},{\tilde  {\rho _{0}}} because of the different expressions for s_{0},g_{0},{\tilde  {s_{0}}},{\tilde  {g_{0}}}. However, the overall governing principle still remains. Therefore, both the original and modified model are subject to the same instability conditions.


The main findings of this investigation are:

  • Populations should be kept below a certain critical density threshold that may be calculated based on known biological parameters, A_{g},a_{g},R_{{s,g}},r_{{s,g}}.
  • Choosing exponential functions for the phase conversion rate expressions yields the same instability condition as for rational phase conversion rate expressions.
  • Instability conditions are governed by the expressions for social interaction potentials Q_{{s,g}}.


This section contains MATLAB code to produce the plots on this page.

close all

%Spatially homogeneous steady states
k1 = 20; %locusts/m^2
% k1 = 65;
k2 = 65; %locusts/m^2
d1 = 0.025; %hr^-2
% d1 = 0.25;
d2 = 0.25; %hr^-2
p0 = 1:1e4; %locusts/m^2
%original model
f1 = d1 ./ (1 + (p0*(1/k1)).^2);
f2 = d2*(p0*(1/k2)).^2 ./ (1 + (p0*(1/k2).^2));

%modified model
f11 = d1 ./ (1 + exp(p0/k1));
f22 = d2*exp(p0/k2).^2 ./ (1 + exp(p0/k2));

% loglog(p0,f1)
% hold on
% loglog(p0,f2,'g')
% legend('f_1','f_2','Location','NorthWest')
% title('Density Conversion Rates - Original Model')
% xlabel('\rho_0 (locusts/m^2)')
% ylabel('f_1, f_2')
% figure
% loglog(p0,f11)
% hold on
% loglog(p0,f22,'g')
% legend({'$\tilde{f_1}$', '$\tilde{f_2}$'},'interpreter', 'latex','Location','NorthWest')
% title('Density Conversion Rates - Modified Model')
% xlabel('\rho_0 (locusts/m^2)')
% ylabel('$\tilde{f_1}$, $\tilde{f_2}$','interpreter', 'latex')
%axis([0 10 0 0.5])

%original model
s0 = (p0 * d1 * k1^2 .*(k2^2 + p0.^2)) ./ ...
    (d1*k1^2*(k2^2+p0.^2) + d2*p0.^2.*(k1^2 + p0.^2));

g0 = (d2 * p0.^3 .* (k1^2 + p0.^2)) ./ ...
    (d1*k1^2*(k2^2+p0.^2) + d2*p0.^2.*(k1^2 + p0.^2));

%modified model
s00 = (p0*d1.*(1 + exp(p0/k2))) ./ ...
    (d1*(1+exp(p0/k2)) + d2*(exp((p0/k1)+(2*p0/k2)) + exp(2*p0/k2)));

g00 = (p0*d2.*exp(2*p0/k2).*(1 + exp(p0/k1))) ./ ...
    (d1*(1+exp(p0/k2)) + d2*(exp((p0/k1)+(2*p0/k2)) + exp(2*p0/k2)));

% figure
% loglog(p0,s00,':k')
% hold on
% loglog(p0,g00,'k--')
% xlabel('\rho_0 (locusts/m^2)')
% ylabel('$\tilde{s_0}$, $\tilde{g_0}$','interpreter', 'latex')
% title('Spatially Homogeneous Steady States, k_1 = k_2 and \delta_1 = \delta_2')
% legend({'$\tilde{s_0}$', '$\tilde{g_0}$'},'interpreter', 'latex','Location','NorthWest')
% axis([0 10 0 10])

hold on
xlabel('\rho_0 (locusts/m^2)')
ylabel('s_0,g_0 (locusts/m^2)')
title('Spatially Homogeneous Steady States, k_1 \neq k_2 and \delta_1 \neq \delta_2')
legend('s_0','g_0', 'Location','NorthWest')

External links


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