
MBW:Two Models of Locust Phase ChangeFrom MathBioContentsExecutive SummaryIn 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 reactiondiffusion 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:
Context/Biological phenomenon under considerationLocusts 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. HistoryFor 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 (integroPDEs). Below, the classic swarm model is described. The findings of this investigation are based on a modified version of this model. Assumptions
Model constructionClassic swarm modeling involves a conserved population density field travelling at a velocity that arises from social interactions:
The contribution of of a small group of individuals at location to the force on the individual at position is given by:
Therefore,
The reader should note that the term is radially symmetric. Let us use some notation to represent this integral:
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
Full modelAssumptions
Model constructionThere are separate density fields for the solitarious and gregarious phases , and these phases sum to the total density:
The reactiondiffusion equations of the model are:
Where the velocities are:
We now require equations for the solitarious/gregarious social interactions and the two densitydependent conversion rate functions, , which represent the rate of gregarious solitarious and the rate of solitarious gregarious, respectively:
Failed to parse (lexing error): Q_g(\mathbf{ xx'}) = R_g e^{\mathbf{xx'}/r_g} – A_g e^{\mathbf{xx'}/a_g}
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, . The modified functions are:
The purpose of introducing exponential functions 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 and are monotonically decreasing and increasing, respectively. Critical distanceFinding the critical distance leads to certain requirements in the model with respect to the interaction amplitudes and interaction length scales, . These calculations are based on the social interaction potentials, . 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 is a decaying exponential function for all choices of and , so it is purely repulsive. However, , being the difference of two exponentials, suggests that there is an equilibrium, where there is no net contribution to the velocity. The equilibrium point is found by solving .
Failed to parse (lexing error): R_ga_g – A_gr_g > 0
Failed to parse (lexing error): A_ga_g^2 – R_gr_g^2 > 0
ParametersThe numerical values of these parameters are based on estimates from biological literature ^{[1]}.
Analysis/interpretationThe analysis of these two models has two parts:
Homogeneous steady statesThe steady states of the model can be obtained by setting the derivatives in time and space to zero and using the expression for .
Thus, the expressions for the equilibrium points are:
For the original model, the full expression is:
For the case where , this is a plot of the solutions: For the case where , this is a plot of the solutions: The corresponding equations for the modified model:
For the case where , this is a plot of the solutions: The following is a zoomed in plot near the origin: For the case where , this is a plot of the solutions: The following is a zoomed in plot near the origin: Linear stabilityThe process of performing a linear stability analysis on the HSS solutions is outlined below:
Consider small perturbationsThe perturbations about the HSS solutions are denoted as :
Accordingly, our expression for changes slightly.
For convenience, refer to this table that shows the analysis of the original model and the modified model side by side:
LinearizeSubstitution and further simplification of the perturbed equations into the reactiondiffusion equations leads to a linearized system. For reference, the reaction diffusion equations are:
Fourier expansionA Fourier expansion is required to further analyze the linear system. We expand the perturbations as an infinite Fourier series. The domain is infinitely large. is the wave number.
Construct matrix of each Fourier mode amplitudeUsing the Fourier expansions in the linearized system reveals a ODEs for each Fourier mode amplitude. The Fourier transforms of the social interaction potentials, are defined as:
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:
Eigenvalues determine instability conditionsThe 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.
Since are always negative, we look to for the instability conditions. In other words, we are looking for the point where . At this point, let us rewrite as mass fractions that sum to 1 (for convenience).
now become:
The point where 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 of the expression
leaves us with:
The expressions for are essentially identical at this point, excepting the terms. Furthermore, the two functions increase monotonically with , 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 . For the first term in the square brackets to be small, the factor of must be small. This implies that 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 and . Therefore, the first two terms of the bracketed quantity are either constant or monotonically increasing. Therefore, the minimum occurs where . Using in the expression for 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 }
DiscussionThe main findings of this investigation are:
CodeThis section contains MATLAB code to produce the plots on this page. clear clc 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]) 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') External linksReferences
