May 21, 2018, Monday

# MBW:Two Models of Locust Phase Change

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

Locust life cycle [3]

## History

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.

### Assumptions

• 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)$

Therefore,

${\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).

### Parameters

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

### Assumptions

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

$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}}}$

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}}}}}$

Plot of ${\tilde {f}}_{{1,2}}(\rho )$
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

### Parameters

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

## Analysis/interpretation

The analysis of these two models has two parts:

• Linear stability

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:

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:

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:

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

The following is a zoomed in plot near the origin:

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:

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

The following is a zoomed in plot near the origin:

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

#### Linearize

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:

$\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.

## Discussion

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}}$.

## Code

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

clear
clc
close all

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

figure
loglog(p0,s0,':k')
hold on
loglog(p0,g0,'k--')
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')