May 23, 2018, Wednesday

# Gypsy moth population dynamics

## Article

Following a paper from J.W. Wilder about gypsy moth population dynamics [1] will be discussed.

## Background

Since gypsy moths are one of the most voracious forest pests discovered, a model of the life cycle might explain the breakouts which are responsible for their amazing destruction. The life cycle of them can be separated in different stages, which can additionally be divided into age distinct cohorts. Every stage has different predators and parasites and the growth of the stages also depends on weather, climate and available Food (foliage). Fig.1 below gives an idea about the relation of all the factors mentioned.

File:Push arrows.jpg
Fig.1: Relation of all influences

Considering all these parameters in a model is nearly impossible. An assumption has to be made for a simple model. The following model is such an assumption. As noticed in the paper the model will fit the empirical data very well.

## Model

As already mentioned the model is a simple - 2 dimensional - model. The two differential equations describe the biomass densities of gypsy moth larvae (G) and foliage (F) (the value and description of each parameter is shown in tab. 1):

Ne is defined as the total number of eggs laid per ha. As shown in (3) it depends on the population density G.

Because of the hatched larvae of the survived eggs and because of the amount of survived larvae the initial condition for the larvae biomass density has to be recalculate after every year, which is shown in (4). A similar dependency of the influences of the foliage density results also in calculating a new initial condition for F every year - shown in (5).

Fcp is representing the growth rate of foliage, where cp is the carrying capacity of foliage of the stand. This capacity is recalculated after every year because of the dependency on the consumed foliage from the larvae the last year. fd and fr give an idea about the decrease of carrying capacity in case of total defoliation (fd) respectively fr calculates the recovery of carrying capacity if there is no more defoliation.

In (13) and (14) the first initial conditions for the simulation start are shown.

Parameter Value Description
Ws 1.0 mass (g) of a pair of large larvae in a sparse population
Wd 0.24132 mass (g) of a pair of large larvae in a dense population
wL 2.0x10-4 mass (g) of a newly hatched larva
sw 0.70 percent of larvae surviving the winter
sd 0.87 percent of larvae surviving dispersal
Ns 425 eggs/egg mass in a sparse population
Nd 272 eggs/egg mass in a dense population
pd 0.8309 percent decrease in carrying capacity
nd 1.0 number of years of defoliation resulting in pd
nr 5.5 number of years required for recovery of stand
Tlife 645 length of larval growth cycle in degree days
fi 43.8 factor by which small populations increase per year
File:Alpha new.png 0.10 assimilation efficiency of larvae
cpmax 2x106 maximum carrying capacity of stand (g per ha)
n 4.01 decay rate of larvae due to severe overpopulation
P 850 predator biomass density (g per ha)
c2 1.51768x10-5 g of larvae consumed per g of predator per degree day
c5 3.33x10-3 search rate of larvae for foliage
c6 3.03687x10-2 relative growth rate of foliage (degree day-1)

Tab. 1: Parameters of the model[2]. The values of fi and cpmax varied during the simulations.

## Results

In the paper the resulting plots were fitting the data very well. To reach this two parameters - fi and cpmax - have been varied in defined ranges (20 <= fi <= 55, 1.5x106 <= cpmax <= 2.5x106). In Fig.1 is a the resulting plot shown, which fits the data. Each peak in the plot stands for an outbreak.

 300px 300px Fig.2: Simulations which fit data with one outbreak [3] Fig.3: Simulations which fit data with two outbreaks [3] 300px 300px Fig.4: Simulations which fit data with two outbreaks [4] Fig.5: Simulations which fit data with two outbreaks [5]

In Fig.6-9 are plots shown, which are the result of the generated matlab function. The curves are not completely similar, because the exact parameters settings of Fig.2-5 were not mentioned in the paper, but they give an idea about the behaviour of the model.

 350px 350px Fig.6: Simulation: one outbreak with fi=30, cpmax=1.5x106 Fig.7: Simulation: two outbreak with fi=45, cpmax=2.0x106 350px 350px Fig.8: Simulation: three outbreaks with fi=37.1429, cpmax=1.5x106 Fig.9: Simulation: four or more outbreaks with fi=38.5714, cpmax=1.5x106

## Matlab Code

### Main code

The code belowing is showing all plots for a varying fi in the range mentioned before.

Note: To run the code the function GFrhs is required which is shown at the bottom of this chapter.

```%% By Stephan Kupsa

global alpha P cpmax cp n sw sd Ns fi Tlife Ws Wd Ns Nd ll ul m wL pd nd nr c1 c2 c3 c4 c5 c6 fr fd cc

alpha = 0.10;
P = 850;
cpmax = 2.5e06; % or 1.5e06 or 2.0e06
cp = cpmax;
n = 4.01;
sw = 0.70;
sd = 0.87;
fi = 43.8;
Tlife = 645;
Ws = 1.0;
Wd = 0.24132;
Ns = 425;
Nd = 272 ;
ll = 1.25e04 * Ws/2; % lower limit of biomass density
ul = 10e06 * Wd/2;   % upper limit of biomass density
m = ul  - ll;
wL = 2.0e-04;
pd = 0.8309;
nd = 1.0;
nr = 5.5;
c1 = n/(Tlife*log10(exp(1)));
c2 = 1.51768e-05;

c5 = 3.33e-03;
c6 = 3.03687e-02;
fr = (1/((1 - pd)^(1/nr)))-1;
fd = (1-((1 - pd)^(1/nd)));

kk=0;% counter1 for the loops
for fi=linspace(20,55,50)
G0 = 2 * wL;
c3 = alpha^(-1);c4 = log(sw*sd*Ns/(2*fi)) / (P*Tlife);
kk=kk+1;
fi_(kk)=fi;
cp = cpmax;
F0 = cp/498.5;
cc = 0; % counter2 for the loops
for t_=1:20

[t,GF]=ode15s(@GFrhs,[0,Tlife],[G0;F0]);

if GF(length(GF),1) <= ll
Ne = GF(length(GF),1) * Ns / Ws; % total # of eggs laid per ha
elseif ll < GF(length(GF),1) <= ul
Ne = GF(length(GF),1) * ((m*Ns-(Ns - Nd)*(GF(length(GF),1) - ll))/(m*Ws-(Ws - Wd)*(GF(length(GF),1) - ll)));
elseif GF(length(GF),1) > ul
Ne = GF(length(GF),1) * Nd / Wd;
end

Fcp=(F0*cp)/(F0+(cp - F0) * exp(-c6 * t_));

if 1 - GF(length(GF),2)/Fcp <= fr
cp = cp*(fr + GF(length(GF),2)/Fcp);

elseif 1 - GF(length(GF),2)/Fcp > fr
cp = cp * (1 - fd * (1 - fr)^(-1)*(1 - fr - GF(length(GF),2)/Fcp));
end

cc = cc + 1;
G0 = wL * sd * sw * Ne;
F0 = cp/498.5;
x(cc)= t_;
y(cc)= Ne;

end
figure(kk);
y_norm = y./max(y);

plot(x,y_norm);
title(num2str(fi));
xlabel('time(years');
ylabel('Normalized Egg Mass Density');

end
```

### Function code

```function GFdot=GFrhs(t,GF)

global alpha P cpmax cp n sw sd Ns fi Tlife Ws Wd Ns Nd ll ul m wL pd nd nr c1 c2 c3 c4 c5 c6 fr fd cc

GFdot(1,1) = -c1 * GF(1)+ alpha * c2 * ((GF(1) * GF(2))/(1 + c2 / c3 * GF(2))) - c4 * ((GF(1) * P)/(1 + c4 / c5 * GF(1)));
GFdot(2,1) = c6 * GF(2) - c6 / cp * GF(2)^2 - c2 *((GF(1) * GF(2))/(1 + c2 / c3 * GF(2)));
```

## References

1. Wilder, J.W., 1999, A predictive model for gypsy moth population dynamics with model validation, Ecological Modelling, 116:165–181,
2. Wilder, J.W., 1999, A predictive model for gypsy moth population dynamics with model validation, Ecological Modelling, 116:179-180
3. Wilder, J.W., 1999, A predictive model for gypsy moth population dynamics with model validation, Ecological Modelling, 116:174
4. Wilder, J.W., 1999, A predictive model for gypsy moth population dynamics with model validation, Ecological Modelling, 116:175
5. Wilder, J.W., 1999, A predictive model for gypsy moth population dynamics with model validation, Ecological Modelling, 116:176