Population dynamics of the forest tent caterpillar (Malacosoma disstria) in a water tupelo (Nyssa aquatica) forest: A simulation model

Ecological Modelling, 39 (1987) 287-305

287

Elsevier Science Publishers B.V., Amsterdam - Printed in The Netherlands

POPULATION

DYNAMICS OF THE FOREST TENT CATERPILLAR

(MALA COSOMA DISSTRIA) IN A WATER TUPELO ( N Y S S A A QUA TICA) F O R E S T : A S I M U L A T I O N M O D E L

MARCEL REJMANEK :, JAMES D. SMITH and RICHARD A. GOYER

Department of Entomology, Louisiana Agricultural Experiment Station, Louisiana State University Agricultural Center, Louisiana State University, Baton Rouge, LA 70803 (U.S.A.) (Accepted 2 October 1986)

ABSTRACT Rejmfinek, M., Smith, J.D. and Goyer, R.A., 1987. Population dynamics of the forest tent caterpillar (Malaeosoma disstria) in a water tupelo (Nyssa aquatica) forest: a simulation model. Ecol. Modelling, 39: 287-305. A population model of the forest tent caterpillar, Malacosoma disstria Hub., which simulates both larval grazing in the forest canopy in one season and population dynamics over several years, is presented. The FTC is an important pest of hardwoods causing very irregular but frequently complete defoliation of water tupelo, Nyssa aquatica L. The FTC population model consists of two submodels: the first one for larval mortality and grazing of growing foliage in the forest canopy, the second for year-to-year population dynamics under different levels of density-dependent pupal mortality due to parasitism. The model prediction of complete defoliation (over 98%) for > 22 egg masses per tree is in agreement with available data. Under conditions of strong density-dependent pupal mortality (rather dry tupelo forests) the modeled population exhibits damped oscillations and approaches equilibrium in a few years. Weak density-dependent pupal mortality, which seems to be typical for permanently flooded forests, results in chaotic or pseudoperiodic model behavior. A method to achieve more regular dynamics through introduction of spatial heterogeneity into the model is illustrated.

INTRODUCTION I n spite of d r a m a t i c effects of s o m e insect species o n d e c i d u o u s forest canopies, there are still o n l y a few a t t e m p t s to use m a t h e m a t i c a l m o d e l s for s i m u l a t i o n of the d e f o l i a t i o n process (e.g., Reichle et al., 1973; V a l e n t i n e et

a Present address: Department of Botany, University of California, Davis, CA 95616, U.S.A. 0304-3800/87/$03.50

© 1987 Elsevier Science Publishers B.V.

288 al., 1976; Valentine and Talerico, 1980; Valentine, 1981, 1983; Rubtsov and Shvytov, 1980; Sollins et al., 1981) or long-term population dynamics of leaf-consuming insects (e.g., Morse and Simmons, 1979; Brown et al., 1983; Rubtsov, 1983). Our site model of the forest tent caterpillar (Malacosoma disstria Hiibner, Lepidoptera: Lasiocampidae) combines both the seasonal modelling of larval grazing in a forest canopy and the long-term, year-to-year modelling of population dynamics. The forest tent caterpillar (FTC) is a native pest of hardwoods throughout most of the United States and Canada. In the South, water tupelo (Nyssa aquatica L.) is the preferred host. Defoliation of this resource causes disruption of the swamp environment and leads, subsequently, to reduced tree growth and, presumably, to tree mortality. Population studies of the FTC have been conducted in the Lake States where trembling aspen (Populus tremuloides Michx.) was the predominant host (Witter et al., 1972). Studies involving water tupelo (e.g., Oliver 1964; Abrahamson and Harper 1973; Stark and Harper, 1982) have been centered on identifying mortality agents with little, or no effort on population modelling. Our study, an extension of the population dynamics study of Smith (1983), is designed to simulate FTC population in permanently flooded water tupelo forest and in areas where standing water occurs only during the springtime. These two hydrologic regimes are typical for water tupelo dominated forests in Louisiana and in some other coastal states. STRUCTURE OF THE MODEL The FTC model consists of a set of differential equations and certain ancillary functions that provide either initial values of the components or values of parameters. Although several facts, assumptions and techniques are adapted from other authors, the major data base is extracted from Smith's (1983) study of the FTC in south Louisiana.

Defoliation process and larval population Defoliation is a central process in the system examined. It proceeds during 29 days of larval development in the canopy of water tupelo. The percentage of defoliation determines the incremental growth loss of trees while the degree of defoliation depends upon the larval population as influenced by biotic and abiotic mortality factors. A differential equation describing the growth of foliage being consumed by larval population, L (t) (m 2 ha-l), is adopted from Reichle et al. (1973) and Valentine (1983): dL(t)

dt

L(t) d M ( t ) M(t) dt

- - -

kl c ( t )

(1)

289

where M(t) is the area of foliage (m 2 ha -1) that would exist in the absence of consumption, C(t) the rate at which foliage is being consumed at time t (m 2 day-a), and k 1 is a constant which represents the area of foliage both dropped and consumed, divided by the area consumed (kl = 1.05 in all our simulations). In agreement with Rubtsov and Shvytov (1980) and Valentine (1983) we suppose a logistic potential growth of foliage, M ( t ) : dM(t)

dt

= 1"1 m ( t ) (1 - m ( t ) / k 2 )

(2)

where r~ is relative growth rate of foliage (m 2 m -2 d a y - a ) , and k 2 maxim u m possible leaf area (m 2 ha-a). On the basis of available data from representative water tupelo forests in south Louisiana we use r I = 0.26 m 2 m -2 day -1, k 2 = 35000 m 2 ha -1 and M(0) = L(0) = 200 m 2 ha -a as typical values. If C(t) = 0 for every t > 0, then L(t) = M ( t ) ; therefore solutions of equations (1) and (2) are identical. The percentage of defoliation at the beginning of pupation is calculated as 100(1 - L(t)/M(t)), t = 29. The consumption rate, C(t), is a product of the mean actual individual consumption, Ale(t) (m 2 d a y - l ) , and the total n u m b e r of larvae per ha, N(t):

C(t)=AIC(t) U(t)

(3)

The mean actual individual consumption, Axe(t), is defined as: AIe(t) = Pie(t) (1 - e x p [ - k 3 L(t)/M(t)]}

(4)

where Pie(t) is potential individual consumption (see Smith et al., 1986) and k 3 is a constant representing searching efficiency. The dependence of the A I e / P I e ratio on the L / M ratio for different k 3 values is illustrated by Fig. 1. The true value of k 3 must be assessed experimentally; k 3 = 8.0 is used in all presented simulations. The A I c / P I e ratio depends, in our opnion, rather on the L / M ratio than on the ratio L / ( t o t a l larval biomass per ha) as suggested, e.g., by Rubtsov and Shvytov (1980) and Valentine (1981). Even if the total larval biomass is low compared to L, the small L / M ratio, say less than 0.1, means that accessibility of foliage is low because its density in the canopy is low. This is correctly reflected by function (4). The dependence of the P I e ( t ) o n the weight of individual larvae, W(t) (g), is described by an empirical S-shaped function based on rearing experiments:

PIC(t)=k4/{1 + k 5 exp[-k 6 W(t)] }

(5)

where k4, k 5 and k 6 are constants estimated by nonlinear regression analysis as 1.93 × 10 -4, 76.2 and 1800, respectively (Fig. 2). (The BMDP

290 LO.

lo/~,//~,.~ .."

0.8-

//

~

/

/ / I"

I/,;

]1,"

0.6.

O

K

0.4-

<

0.2-

I/ / /I /

s

0.0-

o.o

o'.2

o;8

o14

o'.8

alo

I/M Fig. 1. Actual individual consumption/potential individual consumption ratio, A I C / P I C , versus actual area of foliage/maximum area of foliage ratio, L/M, for three different levels of searching efficiency, k 3.

program PAR was used for all nonlinear regression analyses; Dixon 1983). The growth of an average larva is described by the differential equation: dW(t) dPw(t) d-----~ - sup d~

(6)

where SUP =

k 7 R / ( R + k8)

- k9,

R =

AIC(t)/PIC(t)

(7)

and a w ( / ) is potential weight, i.e., mean weight of an individual larva at time t if AIC PIC in a time interval [0, t]. The equation (7) represents the functional response of larvae to a varying energy resource base. Parameters =

0.0020-

0.00015. o

"o

0.000~

vE

tO O.

0.0005.

0.0000-

-4.0

-315

-3'.0

-Ls -2'.0

-,'.5

-I.'O

log 10 W(O)

Fig. 2. Nonlinear regression of the potential individual consumption, PIC, on the weight of individual larvae, IV, based on rearing experiments.

291 tOO, 0.75 0.~,

0.25. n

0.00

tt)

-0.25, -0.~.

oo

dz

o.'4

o76

o'.a

(o

AIC/PIC Fig. 3. Dependence of the parameter suP (equation 6) on actual individual consumption/ potential individual consumption ratio, AIC/PIC.

4.0, k 8 = 2.0, k 9 = 0.33) were chosen to be in agreement with pertinent data in the literature (Gutierrez et al., 1984). The ratio a I c / P I C = 0.18 represents the compensation point (suP = 0; energy captured in food equals the cost of respiration). Fig. 3 shows the dependence of suP on the ratio A I C / P I C . Laboratory experiments are required for better assessment of parameters in equation (7). Dynamics of the individual potential weight is given by empirical function: (k 7 =

dPw(t) dt

( 0, = , r2 PW(t)(1-PW(t)/klo),

0 < t < 3.0 3.0
(8)

where r2 = 0.39 is relative growth rate in g g-1 day-a estimated by nonlinear regression analysis (Fig. 4), and kl0 = 0.083 g is maximum mean larval weight at the end of the 5th instar. Mean larval weights from rearing experiments (Fig. 4) were in close agreement with low density field data (Smith, 1983). The FTC larval population declines from the time of egg-hatch until pupation: dN(t)

d----~ = --MOR

N(t)

(9)

where MOR = 0.24G + 2.16G 2 - 1.44G 3 + 0.04, in which: G = 1 -

(W(t)/PW(t))

2

(10)

Mortality, MOR, then is dependent upon the ratio of mean actual weight to potential weight of larvae at time t. This ratio integrates the history of starvation and therefore is less mechanistic than the ratio (total larval

292

0.08-

0.07. ¢~ II I.d

~=

0.060.050.04.

0.o5Iz LU I.0 n

0.020.01-

O.O0-

TIME ( D A Y S )

Fig. 4. Average potential weight of larvae, PW, versus time. The data from laboratory experiments were fitted by equation (8).

potential consumption)/L which is usually used for the same purpose. Function (10) is depicted in Fig. 5. The form and parameters of the function (10) were chosen to be in best possible agreement with our data (0.04 is the mean field mortality for densities less than 1 million of hatched larvae per ha) and with available information about survival under food deprivation of Lepidoptera larvae (Hodson, 1941; Beckwith, 1983). More laboratory and field work is required to improve the starvation-mortality relationship. Equations (1) through (10) are solved simultaneously in the course of time interval 0-29 days starting with synchronized water tupelo budbreak and hatching of the FTC each year.

[.0" 0.8"

0.6" n0

0.4-

0.2-

0.0-

oo

o'2

o'4

o.'6

o'8

,Jo

1 - W/PW

Fig. 5. Assumed dependence of the mortality parameter, MOR, in equation (9) on the ratio mean actual larval weight/potential larval weight, W/Pw.

293

Year-to-year population dynamics In the year-to-year population dynamics, final larval density (initial pupal density), pupal mortality and presumably density-dependent fertility play the most important roles. The pupal mortality is given by Hassell's (1975) density-dependent function: ADU = PUP(1 + a a

(11)

PUP) -~

where ADU is number of adults per ha, PUP is initial number of pupae per ha (PEr' = N(t), t = 29) and a a and fl are empirical constants estimated on the basis of available data (Fig. 6). Equation (11) can be rewritten in the form 1og(r'UP/ADU) = fl log(1 + a a PUP) which is used, e.g., by Hassell et al. (1976). Parameter fl is the slope of the relationship attained at high population densities and a 1 is related to the point of inflection of the curve. Sarcophaga houghi Aldrich, a sarcophagid parasitoid, is the most important agent causing mortality of pupae. Our data (Fig. 6) suggest two different intensities of parasitoid control: weak on the localities 'Alligator' and 'Verret' (A), and rather strong on the locality 'Sorrento' (B). The situation A represents permanently flooded forests; in B standing water occurs only during the spring time. Even if the role of flooding in the regulation of S. houghi was questioned (Stark and Harper, 1981), our data support rather the conclusions of earlier authors who attributed the lack of parasite control of FTC to flooding which drowns overwintering puparia of S. houghi and other parasitoids (Hodson, 1941; Abrahamson and Harper, 1973, Batzer and Morris, 1978).

3.0 J r~ P

2.5-

_o (.9 0 /

)_J I--

o J

O_

2.0-

B ,'" f

1.5L00.50.0-

i

1.8

2.4

3.0

3.6

4.2

4.8

5.4

6.0

PUPAE PER HECTARE, LOGI0 (PUP)

Fig. 6. Density-dependent pupal mortality data fitted by equation (11) in habitats A and B. PuP is initial number of pupua per ha, ADU is number of adult per ha.

294

Assuming a sex ratio = 0.5, the number of eggs masses per ha, Mh, is simply: M h = 0.5 ADU

(12)

Mean number of egg masses per tree, MT, then depends upon the n u m b e r of trees (dbh > 3 in = 7.6 cm) per ha, T:

MT = M h / T

(13)

In our plots there were 667 (Sorrento), 666 (Verret) and 882 (Alligator) trees of the diameter 3-32 in (7.6-81 cm) per ha. We use, therefore, T = 700 for all simulations. Fecundity, or mean number of eggs per egg mass, EM, is considered either constant, namely E M = 346.5 eggs per mass, or dependent on the initial weight of pupa according to the equation: EM=

370(W(I)/Pw(I)),

t = 29

(14)

The number of eggs per ha, Eh, is given as: (15)

E h = EMM h

and survival of eggs is assumed to be constant:

No=N(t)=O.778E h,

t=0

(16)

Equation (16) defines initial larval density at the time of budbreak. Subsequent leaf consumption and larval population dynamics are defined by equations (1) to (10). The program is written in SAS language using DATA and PROGRAM statements. Procedures MODEL and SIMNLIN f r o m SAS/ETS Library were used for numerical solution of differential equations. RESULTS OF SIMULATIONS

Defoliation process and larval population dynamics The simulation of the foliage growth and defoliation in a single season is shown in Fig. 7 for three different initial densities of egg masses per tree, M T, provided 346.5 eggs per mass. F r o m our simulations it follows that M T less than 5 does not cause any substantial damage (less than 26% defoliation). The 50% defoliation corresponds to M T = 10.1. The M T greater than 15 causes defoliation greater than 70%. Complete defoliation (greater than 98%) starts with M T greater than 22. Results of simulations with parameters specified above are in agreement with available M T data and estimates of defoliation (Smith, 1983). Corresponding changes in larval density are depicted in Fig. 8. The mortality due to the starvation starts to be important for M T higher than 15.

295 40000~

A

30000 5

RO000-

i.iJ rr I0000-

15

Ld ¢D

8

o

I.L

;

,b

,;

zb

£5

3'o

TIME (DAY) Fig. 7. Simulation of the foliage growth for four different initial densities of egg masses per tree ( M T = 0, 5, 15, 20), provided 346.5 eggs per mass.

The relationships between initial density of 1st instar larvae, No, and resulting pupal density PuP, is shown in Fig. 9 for three different relative growth rates of foliage (r I = 0.20, 0.26, 0.32). Results of simulations based on equation (3) are indicated as crosses and are fitted by function (17). Empirical data from Smith (1983) are indicated as closed circles. Simulations with r I = 0.26 are in agreement with available data.

Year-to-year population dynamics with constant fecundity Long-term population dynamics of FTC was simulated for two different combinations of parameters a 1 and fl in equation (11), i.e., for two different intensities of density-dependent pupal mortality. Values a 1 = 1.9 × 10 - 4 ,

6.50

.oo1 1 1 i\ _~

5.254,

o

I

I

,I

I

,

~

lb

l~

20

z~

' 30

TIME (DAYS) Fig. 8. Decline of larval densities for three different intial m e a n densities of egg masses per tree ( M T = 5, 15, 20), provided 346.5 eggs per mass.

"~l!punooj lUl~lSUOO Ioj V lgl!ql~q oql u.t uo!ll~.loJo p j o p u e 'I/4r 'ooal aod sossem ~'~o j o s~oqtunu u~otu j o so.n.ueu~p m~ol-~uo I go uo.tl~Inm!s "0[ "~!fl

(l:lV3~,) 31~I11 I'll Og i

g~ .

O~ g£ , i

0.£ .g~ . O~. gl 0 . .OI g ============================================== -0

::0 ZO

:::::::::::::::::::::::::::::::::::::::::::::::: -0~

"Or,

"09

-og

rq

001

s:tolotur.lrd jo oogds oql u! [opom Oql jo s!s~lrur oql ol sn po[ llnso~ luo~ojj!p ~[~u.p~. lS s.tq,.L "uo.tlg!lojop %6"E[ pu~ ooal ~tod sossgtu g~o 9"E lg tun!.~q.q -mbo oql ol qogo~tddg faol~ii!oso u! sllnsoi fl l~l!qgq u! ,q!imaotu lgdnd luopuodop-,(l!suop ~o~uo~tls 'fa~.tluoo oql u o "uo!lglntu!s oql jo un.~ .to~uoI g ,~q pol~:tlsuotuop oq ugo s~ o!po~.~odg s~. so.ttugu£p oql 'logJ uI ",solo~3, ago~-~ luo2gddg sh~ol[~ V lm!qgq u! gl~.lglaotu p~dnd luopuodop-£1t.suop ~IeOA~ "([[ pug 0[ "sg~.d) luo.tojj!p ,([olo[dtuo3 'suo.~l!puoo jo los o~1 osoql ~opun 'oag so.rtueuKp uo!lglndod polgintu!S "(9 "~!3 oos) ~Ioa!loodso.t '~ pue V suo!l -rnl!s o.~oIoap,( q .~oj polrtugso o~o~ 8['~ = g '~_0~ x 9"[ = ~v p u t 91~'0 = 0t

"SOlOIt.oposolo sg polro!pu! olr i~li~p I~oy!dtuo olqg[.reAV "(LD uo.tlounj ~q poilu ose pug sosso~o se poleo!pu! o~g (6) uo!lenbo uo posgq uo!lrlntu!s jo sllnso~ "(EC0 '9E'0 '0I'0 = L0 o~g.qoj jo s o l ~ q l ~ o ~ o^!1~IOl luo~ojj.tp oo~ql Ioj 'oN 'ogmg I aelSUl, ls[ jo aoqtunu ll~!Lmt, oql uo 'aria 'gq aod ogdnd jo aoqtunu oql jo oouopuodo(I "6 "i!fl Dq 83d 3~'A~IV-I ~VJ.SNI J,Sl 40 ~381/~N

Z I

00000001 0000008 0000009 000000~ 000000~ :::::::::::::::::::::::::::::::::::::

ZE 0

'

i

0 •

~,

' 0

00000£

0

000009

C

oy L

9Z 0

0~0

F ~

0

0

0

0

8

~OO

"11

.

1

96E

297

40 an.OZ 5

I0

15 20

25 30

35 40 45

50

TIME (YEAR)

Fig. 11. Simulation of long-term dynamics of mean numbers of egg masses per tree, M r , and of defoliation in the habitat B for constant fecundity.

which determine the intensity of density-dependent relationships involved in the model. Because the original model is too cumbersome for extensive analysis of its behavior in the parameter space, we substitute the defoliation-larval submodel by the function: PuP = cN0(1 +

(azNo)b) -1

(17)

where N o is initial density of 1st instar larvae and c, a2 and b are parameters (see Fig. 9). The density-dependent function (17) was suggested by M a y n a r d Smith and Slatkin (1973). Bellows (1981) compared all suggested density-dependent functions and concluded that function (17) provided the best fit to available data sets. Parameters c (density-independent mortality) and b (severity of the density dependence) are independent on r 1 (specific growth rate of foliage). The value of c follows from the equation (9) and (10): c = 0.313. The value of b was estimated by nonlinear regression analysis: b = 19.8. Only the scaling parameter a 2 depends on r 1. Results of simulations for 20 different r I values in the interval [0, 1.5] were fitted by function: a 2 = 1.19 × 10-8/(2.08 × 10 -3 + rl2"248)

(18)

A simplified version of the long-term population model with constant fecundity than can be written as: ADUt+ 1 = P U P ( 1 + a 1 P U P ) - ~ )

l /

PUP=cN°(1 + (a2N°)b)-I N o = ~ A D Ut

(19)

298

-(~ ADU = PUP (1"oI PUP) .

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

PUP

PUP= CNo(l* (O2No)b)-~

.

l ~-"

b =19.8

k. ADU

==.I===gn

. . . .

No=~ADu $I

= . = . . = = = l l . = = = = = l = l = = = l = . . . =

. . . .

= = . = = l l = = = . =

Fig. 12. To compose the system of difference equations (19), follow the arrows as shown. Two different values of fl represent habitats A and B. The starting point = initial surviving FTC egg density in eggs per ha. After reaching the resultant no. of pupae, stand the graph on its side so that the ' p u p axis' is horizontal. Then, the resulting no. of adults (ADU) depends on the parameter fl used (habitat). The resulting ADU no. is multiplied by )~ to calculate no. eggs surviving.

The parameter )~ involves equations (12), (15) and (16) and its value for constant fecundity, therefore, is 134.8. A system of difference equations (19) is illustrated by Fig. 12 where two different values of fl are used: fig and fib for habitats A and B, respectively. Following the arrows, as shown in Fig. 12, helps to understand qualitative differences in population dynamics in habitats A (Fig. 10) and B (Fig. 11). While population densities in the habitat B are very soon forced to the equilibrium, dynamics in the habitat A is much more complicated (try to follow the process as indicated by arrows). It is clear that in the dynamics of the model defined by (19), the parameters fl and b play a key role. As we may expect somewhat lower values of the parameter b in reality (due to uneven distribution of egg masses per tree) and we have only very preliminary notions about possible values of the parameter fl, the analysis of the behavior of the model (19) in the b-fl space is highly desirable. The results of the numerical analysis are presented in Fig. 13. For this analysis parameters a 1 --1.7 × 10 - 4 (mean value for habitats A and B) and a 2 = 2.35 X 10 -7 (corresponding to r 1 = 0.26) were used. There is very close correspondence between the qualitative behavior of models defined by (1) through (16) and by (19). Slight differences in the quantitative outputs (values of maxima, of equilibria, etc.) are caused by residuals in the regression (17) and by differences between original values of al and mean value of al used in (19).

299

Higher order cycles

2 - p o i n t cycles

Damped oscillations

AB

Monotonic damping

0~hoosDamped oscilations ......

0

5

I0

AA

15

20

25

30

35

b

Fig. 13. Regimes of dynamical behavior of the system (19) as functions of parameters b and B. The shaded area indicates combinations of parameter values producing more than 50% of years with 45% or higher incremental growth loss of trees.

There are various domains of dynamic behavior indicated in Fig. 13. Stability analysis of single functions (11) and (17) have already been done (Hassel et al., 1976; Bellows 1981; see also May, 1983). Complicated behavior of (19), as expressed in Fig. 13, follows from the properties of functions (11) and (17) and for X = 134.8 and fl = 0 or b = 0 is in agreement with published analyses of those functions. Habitats A and B are situated in very different stability domains of b-fl space. From Fig. 13 follows that only slight changes of the fl (down for B, up for A) can change the dynamic behavior dramatically. Bending up of the horizontal fines separating different regions for b < 0.5 in Fig. 13 is produced by function (17). For b = 0, pup = cNo/2, which implies X/2 instead of X in (19). Decrease in ?t changes the boundaries between different regions towards higher values of fl (Hassel et al., 1976, May, 1983). Realism of Fig. 13 for very low b is therefore questionable; however, we do not suppose the parameter b could be less than 1 in any FTC system. The shaded area in Fig. 13 indicates combinations of parameter values producing more than 50% of years with 45% or higher incremental growth

300

loss of trees. It is based on the assumption that two subsequent years with more than 50% defoliation result in 45% incremental growth loss and that defoliation over 50% is caused by densities higher than 10.1 egg masses per tree ( E M = 346.5).

Year-to-year population dynamics with density-dependent fecundity Fecundity, expressed as the mean n u m b e r of eggs per egg mass, EM, is now considered to be negatively dependent on the initial weight of pupae according to the function (14). A n example of simulated population dynamics in habitat A is shown in Fig. 14. Population fluctuations are even more irregular than in the run with constant fecundity (Fig. 10). Also, some values of E M are so low ( < 0.0001) that probably no field population could withstand such decline without becoming extinct. Simulation of population dynamics in habitat B with density-dependent E M produces basically the same result as that in Fig. 11. Only the equilibrium value is slightly higher due to maximum mean fertility (370 eggs per mass) used in equation (14). The relationship between initial density of 1st instar larvae, No, and fecundity, E M, is shown in Fig. 15 (r 1 = 0.26). Results of simulation (triangles) are fitted by the function: EM = g ( 1 + a 3 exp(b3N0)) -1

(20)

where g = 370, a 3 = 1.46 × 10 -5 and b 3 = 1.82 × 10 -5. The only available empirical data from Smith (1983) are indicated as closed circles. So far, the density-dependent reduction in fecundity is not supported by data from water tupelo forests. Also the results of Witter et al. (1975) are not

*

I001 •

bJ L~

i

40-

J

___J

o-

0 ~.h~

I

5

I0 '

.IT.

15 20 25 30 35 40 45 ,

,

,

.

50

,

TIME (YEAR)

Fig. 14. Simulation of long-term dynamics of mean numbers of egg masses per tree, MT, and defoliation in the habitat A for density-dependent fecundity.

301 400.

300

2OO w

,oo. (.3

o.

o

15obooo 3060000 45&~ooo 6o6oooo

NUMBER OF Isf INSTAR LARVAE PER ha

Fig. 15. Fecundity, EM, as a function of the initial density of 1st instar larvae, N0. Results of simulations based on equations (1) through (16) are indicated as asterisks and are fitted by the function (20). Closed circles represent the only available empirical data.

Higher order cycles

2 - point cycles

Damped oscillations

&B Monotonic damping -,~,,,,~,~,,,ffrff ................. ~

i d oscillations X~2- point cycles

' ..... i ..................... 41, A ^.

Tronsffion region

~naoe

0

0

5

I0

15j, 20

25

30

35

U

Fig. 16. Regimes of dynamical behavior of the system (19)-(21) as functions of parameters b and ft. The shaded area indicates combinations of parameter values producing more than 50% of years with 45% or higher incremental growth loss of trees.

302 conclusive in this respect despite a reduction of fecundity to 50% for densities over 20 egg masses per tree reported by Hodson (1941, 1977). For the analysis of the model (19) with density-dependent fecundity, the parameter X was replaced by function: ~. = zg(1 + a 3 exp(b3No)} -1

(21)

where z is the product of sex ratio and egg mortality: 0.5 × 0.778 = 0.389. The behavior of this model is the b-/3 parameter space is shown in Fig. 16. Qualitative behavior of the model (19)-(21) is similar to that of (19), especially for higher values of b. DISCUSSION On the basis of preliminary estimates of several parameters and site models consisting of equations (1) through (16), (19) and (19)-(21), we distinguished two very different kinds of behavior in long-term FTC population dynamics. First (A), where density-dependent pupal mortality is very weak, and strongly nonlinear, density-dependent larval mortality produces chaotic population dynamics. Second (B), where moderate density-dependent pupal mortality does not allow serious density-dependent larval mortality due to starvation and population density soon approaches stable equilibrium (see Figs. 12, 13 and 16). Neither of these two cases corresponds precisely with the available empirical data. However, the model predictions of egg density were verified in 1984 and 1985 by extensive sampling in both habitats A and B. There are greater population fluctuations in habitat A (permanent flooding) comparing with habitat B, but never with such low egg and larval densities following population outbreaks as those shown in Figs. 10 and 14. On the other hand, available data from habitat B (spring flooding) indicate higher population constancy close to the predicted equilibrium value of 2.6 egg masses per tree (Smith, 1983, and unpublished data), but never absolute constancy (cf. Fig. 11). Population fluctuations in habitat B must be caused by factors which are not involved in the model (e.g., climatic fluctuations, massive immigration, or bird predation). Unrealistic minima in the simulated population density in habitat A are in nature undoubtedly buffered by immigration of adults. Our model, therefore, represents a hypothetical situation in a closed system. Any large forest area dominated by water tupelo can be considered as consisting of several patches with different FTC densities at any one time. The number of adults migrating between patches usually is proportional to the population densities in particular patches. Thus the extension of the

303

I.u o3

60

0=

40-~

!

~

20

,

03

z

0

,,,

I'

,,

f 5

I0

15

20

25

30

Z I,LI

TIME (YEAR) Fig. 17. An example of long-term dynamics of the system (22) with n = 3, h,i = 0.998 and h,j = 0.001 for all j. The temporal changes of mean numbers of egg masses per tree in three patches with mutual exchange of adults are indicated by three different lines.

models (19) and (19)-(21) then can be written as: ADUt+I. i + P U P i ( 1 + a 1 PUP/) -/~ (i = 1 , 2 . . . . . n)

pup i = Noi =

cNoi(1+

(azNo, lb) - '

~ ~jh~j ADUt, y

(22)

j=l

= [

134.8

~tj ~zg(1 + 03 exp(b3N0j)} -1

(constant fecundity) (density-dependent fecundity)

where hii = p < 1, E7¢ i hjj = (1 - p ) and n is the total n u m b e r of patches. Coefficient h,, indicates the proportion of adults which remains in a patch i and coefficient hij is the proportion of adults migrating from patch j to patch i. One simulation example using model (22) is shown in Fig. 17 for density-dependent fecundity, n = 3, hig = 0.998 and hij = 0.001 for all j. The resulting dynamics of (22) depends u p o n the initial densities, and Fig. 17 respresents only one of an infinite n u m b e r of possible outcomes. The general feature of the simulations with sufficiently different initial densities and similar h~d coefficients is that population m i n i m a are reasonable (usually more than 0.1 egg masses per tree). Also regular cycles are c o m m o n and frequently established earlier than in Fig. 17. A comparison of Figs. 14 and 17 provides a good idea what introduction of spatial heterogeneity into the model can bring about. Creation of 'order' from 'disorder' is one possible

304 o u t c o m e . Systematic analysis of the b e h a v i o r of the m o d e l (22) is in progress. N o m o d e l can substitute basic ecological research. O u r s i m u l a t i o n m o d e l s were c o n s t r u c t e d to organize existing d a t a a n d i n d i c a t e f u t u r e r e s e a r c h priorities. F u r t h e r research appears to be p a r t i c u l a r l y n e e d e d o n larval searching efficiency, larval s t a r v a t i o n - g r o w t h a n d larval s t a r v a t i o n - m o r t a l ity relationships. Coefficients /3 a n d b which r e p r e s e n t severity of d e n s i t y d e p e n d e n c e in e q u a t i o n s (11) a n d (17), respectively, a n d d i s p e r s i o n coefficients hij in (22) are of u t m o s t interest. ACKNOWLEDGEMENTS W e t h a n k A l a n A. B e r r y m a n a n d James P. G e a g h a n for discussion a n d criticism d u r i n g the p r e p a r a t i o n o f the m a n u s c r i p t . REFERENCES Abrahamson, L.P. and Harper, J.D., 1973. Microbial insecticides control forest tent caterpillar in southwestern Alabama. USDA For. Serv. Res. Note SO-157, 3 pp. Batzer, H.O. and Morris, R.C., 1978. Forest tent caterpillar. USDA For. Serv. Insect Dis. Leafl. 9, 8 pp. Beckwith, R.C., 1983. The effect of temperature and food deprivation on survival of first-instar douglas-fir tussock moth Orgyia pseudotsugata, Lepidoptera: Lymantriidae. Can. Entomol., 115: 663-666. Bellows, T.S., 1981. The descriptive properties of some models for density dependence. J. Anim. Ecol., 50: 139-156. Brown, M.W., Williams, F.M. and Cameron, E.A., 1983. Simulations on the role of the egg parasite, Ooencyrtus kuvanae (Howard), in the population dynamics of the gypsy moth. Ecol. Modelling, 18: 253-268. Dixon, W.J., 1983. BMDP Statistical Software. University of California Press, Berkeley, CA. Gutierrez, A.P., Baumgartner, A.P. and Summers, C.G., 1984. Multitrophic models of predator-prey energetics. Can. Entomol., 116: 923-963. Hassel, M.P., 1975. Density-dependence in single-species populations. J. Anim. Ecol., 44: 283-295. Hassel, M.P., Lawton, J.H. and May, R.M., 1976. Patterns of dynamical behavior in single-species populations. J. Anim. Ecol., 45: 471-486. Hodson, A.C., 1941. An ecological study of the forest tent caterpillar Malacosoma disstria Hbn., in northern Minnesota. Univ. Minnesota Agric. Exp. Stn. Tech. Bull., 148: 1-55. Hodson, A.C., 1977. Some aspects of forest tent caterpillar population dynamics. Univ. Minn. Agric. Exp. Stn. Tech. Bull., 310: 5-16. May, R.M., 1983. Nonlinear problems in ecology and resource management. In: G. Iooss, R.H.G. Helleman, and R. Stora (Editors), Chaotic Behaviour of Deterministic Systems. North-Holland, Amsterdam, pp. 514-563. Maynard Smith, J. and Slatkin, M., 1973. The stability of a predator-prey system. Ecology, 54: 384-391. Morse, J.G. and Simmons, G.A., 1983. Simulation model of gypsy moth introduced into Michigan forests. Environ. Entomol., 8: 293-299.

305 Oliver, A.D., 1964. Control studies of the forest tent caterpillar, Malacosoma disstria, in Louisiana. J. Econ. Entomol., 57: 157-160. Reichle, D.E., Goldstein, R.A., Van Hook, R.I. and Dodson, G.J., 1973. Analysis of insect consumption in a forest canopy. Ecology, 54: 1076-1084. Rubtsov, V.V., 1983. Mathematical model for development of leaf-eating insects oak leafroller taken as an example. Ecol. Modelling, 18: 269-289. Rubtsov, V.V. and Shvytov, I.A., 1980. Model of the dynamics of the density of forest leaf-eating insects. Ecol. Modelling, 8: 39-47. Smith, J.D., 1983. Analysis of forest tent caterpillar, Malacosoma disstria Hbn. Lepidoptera: Lasiocampidae, populations in Louisiana, Ph.D. Diss., Louisiana State University, Baton Rouge, LA, 126 pp. Smith, J.D., Goyer, R.A. and Woodring, J.P., 1986. Instar determination and growth and feeding indices of the forest tent caterpillar, Malacosoma disstria (Lepidoptera: Lasiocampidae), reared on tupelo gum, Nyssa aquatica L. Ann. Entomol. Soc. Am., 79: 304-307. Sollins, P., Goldstein, R.A., Mankin, J.B., Murphy, C.E. and Swartzman, G.L., 1981. Analysis of forest growth and water balance using complex ecosystem models. In: D.E. Reichle (Editor), Dynamic Properties of Forest Ecosystems. Cambridge University Press, Cambridge, pp. 537-566. Stark, E.J. and Harper, J.D., 1982. Pupal mortality in forest tent caterpillar Lepidoptera: Lasiocampidae: causes and impact on populations in southern Alabama. Environ Entomol., 11: 1071-1077. Valentine, H.T., 1981. A model of oak forest growth under gypsy moth influence. In: C.C. Doane and M.L. McManus (Editors), The Gypsy Moth: Research Toward Integrated Pest Management. USDA For. Serv. Tech. Bull. 1504, pp. 50-61. Valentine, H.T., 1983. Budbreak and leaf growth functions for modelling herbivory in some gypsy moth hosts. For. Sci., 29: 607-617. Valentine, H.T. and Talerico, R.L., 1980. Gypsy moth larval growth and consumption on red oak. For. Sci., 26: 599-605. Valentine, H.T., Newton, C.M. and Talerico, R.L., 1976. Compatible systems and decision models for pest management. Environ. Entomol., 5: 891-900. Witter, J.A., Kulman, H.M. and Hodson, A.C., 1972. Life tables for the forest tent caterpillar. Ann. Entomol. Soc. Am., 65: 25-31. Witter, J.A., Mattson, W.J. and Kulman, H.M., 1975. Numerical analysis of a forest tent caterpillar Lepidoptera: Lasiocampidae. Outbreak in northern Minnesota. Can. Entomol., 107: 837-854.