
MBW:Gypsy moth population dynamicsFrom MathBioContentsGypsy moth population dynamicsArticleFollowing a paper from J.W. Wilder about gypsy moth population dynamics ^{[1]} will be discussed. BackgroundSince 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. 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. ModelAs 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): N_{e} 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). F_{cp} is representing the growth rate of foliage, where c_{p} 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. f_{d} and f_{r} give an idea about the decrease of carrying capacity in case of total defoliation (f_{d}) respectively f_{r} calculates the recovery of carrying capacity if there is no more defoliation.
Tab. 1: Parameters of the model^{[2]}. The values of f_{i} and c_{p}^{max} varied during the simulations. ResultsIn the paper the resulting plots were fitting the data very well. To reach this two parameters  f_{i} and c_{p}^{max}  have been varied in defined ranges (20 <= f_{i} <= 55, 1.5x10^{6} <= c_{p}^{max} <= 2.5x10^{6}). In Fig.1 is a the resulting plot shown, which fits the data. Each peak in the plot stands for an outbreak.
Matlab CodeMain codeThe code belowing is showing all plots for a varying f_{i} 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.0e04; pd = 0.8309; nd = 1.0; nr = 5.5; c1 = n/(Tlife*log10(exp(1))); c2 = 1.51768e05; c5 = 3.33e03; c6 = 3.03687e02; 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 codefunction 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
