Forest-Pest Interaction Dynamics: The Simplest Mathematical Models

W O R K I N G P A P E R

FOREST-PEST INTERACTION DYNAMICS:

THE SIMPLEST MATHEMATICAL MODELS

M.

Ya. Antonovekg R.A. Fleming

Yu. A . Kuznetsov W.C. Clark

October 1988 WP-88-092

I n t e r n a t i o n a l I n s t i t u t e for Appl~ed Systems Analysis

FOREST-PEST INTERACTION DYNAMICS:

THE SIMPLEST MATHEMATICAL MODELS

M .

Ya. Antonovsky R.A. Fleming

Yu.A. Kuznetsov W.C. Clark

October 1988 W P-88-092

Working Papers are interim reports on work of the International Institute for Applied Systems Analysis and have received only limited review. Views or opinions expressed herein do not necessarily represent those of the Institute or of its National Member Organizations.

INTERNATIONAL INSTITUTE FOR APPLIED SYSTEMS ANALYSIS A-2361 Laxenburg, Austria

FOREWORD

Some of t h e m o s t exciting c u r r e n t work in t h e environmental sciences involves unprecedentedly close interplay among field observations, r e a l i s t i c but complex simulation models, and simplified but analytically t r a c t a b l e versions of a f e w basic equations. IIASA's Environment Program h a s developed such p a r a l l e l and comple- mentary a p p r o a c h e s in i t s analysis of t h e impact of environmental change o n t h e world's f o r e s t systems.

Two previous p a p e r s (WP-87-70 and WP-87-71) h a v e demonstrated t h e pro- g r e s s t h a t h a s been made. In t h i s new work, t h e conceptual ideas and experimental r e s u l t s contained in t h o s e p a p e r s have been fused together. In p a r t i c u l a r , a sim- ple model of multiple-aged f o r e s t s , t h e i r p r e d a t o r s and t h e i r abiotic environment h a s been developed and successfully t e s t e d with d a t a on budworm populations in North American e a s t e r n s p r u c e forests.

R.E. Munn.

Leader, Environment Program

ABSTRACT

This p a p e r is devoted t o t h e investigation of t h e simplest mathematical models of non-even-aged f o r e s t s affected by insect pests. Two extremely simple situations a r e considered: 1) t h e pest feeds only on young t r e e s ; 2) t h e pest feeds only on old trees. The parameter values of t h e second model are estimated f o r t h e case of bal- s a m f i r f o r e s t s and t h e e a s t e r n s p r u c e budworm. It i s shown t h a t a n invasion of a small number of pests into a steady-state f o r e s t ecosystem could r e s u l t in intensive oscillations of i t s a g e s t r u c t u r e . Possible implications of environmental changes on forest ecosystems are also considered.

SOFIWARE SUPPORT

Software i s available t o allow interactive exploration of t h e m o d e l s described in t h i s p a p e r . The software consists of plotting routines and m o d e l s of t h e systems described h e r e . I t can b e r u n on a n IBM-PC/AT with t h e Enhanced Graphics Display Adapter and 256K graphics memory.

For f u r t h e r information o r copies of t h e software, contact t h e Environment Program, International Institute f o r Applied Systems Analysis, A-2361 Laxenburg, Austria.

-

^vii

-

FUREST-Pm

INTERACTION

DYNAMCS;

THE

SIMP- MA-TICAL

MODELS.

M.Ya. A n t o n o v s k y , R.A. F l e m i n g * , Yu.A. K u z n e t s o v * * and W.C. Clark***

1. Introduction

The influence of i n s e c t p e s t s o n t h e a g e s t r u c t u r e dynamics of f o r e s t systems h a s not been extensively s t u d i e d in mathematical ecology.

S e v e r a l p a p e r s (e-g. Antonovsky a n d Korzukhin, 1983; Korzukhin, 1980) h a v e b e e n devoted to modelling t h e a g e s t r u c t u r e dynamics of a forest n o t a f f e c t e d by pests. Dynamical p r o p e r t i e s of insect-forest systems u n d e r t h e assumption of a g e a n d s p e c i e s homogeneity c a n b e d e r i v e d f r o m t h e t h e o r e t i c a l works o n p r e d a t o r - p r e y system dynamics (May, 1981; Bazykin, 1985). In t h e p r e s e n t p a p e r we a t t e m p t to combine t h e s e t w o a p p r o a c h e s to i n v e s t i g a t e t h e simplest models of non-even- a g e f o r e s t s a f f e c t e d by i n s e c t p e s t s . This p a p e r i s b a s e d upon IIASA WP-87-70 (An- tonovsky e t a l . , 1987); a n d WP-87-71 (Fleming et al., 1987).

The model f r o m Antonovsky a n d Korzukhin (1983) i s a simple model of a g e s t r u c t u r e dynamics of a one-species system. I t d e s c r i b e s t h e time evolution of only t w o a g e c l a s s e s ("young" a n d "old" t r e e s ) . The model h a s t h e following form:

* Forest Pest Management Institute, Canadian Forestry Service, Sault S t e Marie, Ontario, Canada.

** Research Computing Centre. Academy of Sciences of the USSR, Puscheno, USSR.

* 8 8 Science and Public Policy Program, J.F. Kennedy School of Government, Barvard University,

Cambridge, Mass., USA.

where z a n d y are densities of "young" a n d "old" t r e e s , p i s f e r t i l i t y of t h e s p e c i e s , h a n d are d e a t h a n d aging rates. The function 7 ( y ) r e p r e s e n t s a depen- d e n c e of "young" tree mortality o n t h e density of "old" trees. Following Antonovsky a n d K o n u k h i n (1983) w e s u p p o s e t h a t t h e r e e x i s t s some optimal value of "old" tree density u n d e r which t h e development of "young" trees g o e s o n most successfully. In t h i s case i t i s possible to c h o o s e 7 ( y )

=

a ( y

-

^{b ) 2}

+

^c (Figure 1).

Model (A.0) s e r v e s as t h e basis f o r o u r analysis. Let u s t h e r e f o r e r e c a l l i t s p r o p e r t i e s . By s e t t i n g s

=

f + c , scaling v a r i a b l e s ( z , y ), p a r a m e t e r s ( a , b ,c ,p,f , h ,s ) a n d t h e time, system (A.0) c a n b e t r a n s f o r m e d into "dimensionless"

form:

where w e h a v e p r e s e r v e d t h e old notations.

The p a r a m e t r i c p o r t r a i t of system (0.1) on t h e (p,h)-plane f o r a f i x e d s value is shown in Figure 2. Relevant p h a s e p o r t r a i t s are a l s o p r e s e n t e d t h e r e .

Thus, if p a r a m e t e r s ( p , h ) belong t o r e g i o n 2, system (0.1) a p p r o a c h e s a sta- t i o n a r y state with c o n s t a n t a g e c l a s s densities (equilibrium E 2 ) from a l l initial con- ditions. In r e g i o n 1 between lines D l a n d D 2 t h e system d e m o n s t r a t e s a low density t h r e s h o l d : a sufficient d e c r e a s e of e a c h a g e c l a s s l e a d s t o d e g e n e r a t i o n of t h e system (equilibrium E o ) . The boundary of initial densities t h a t r e s u l t in t h e d e g r a - dation i s formed by s e p a r a t r i c e s of saddle El. Finally, in r e g i o n 0 t h e s t a t i o n a r y e x i s t e n c e of t h e system becomes impossible.

Let u s now i n t r o d u c e a n i n s e c t p e s t i n t o model (A.0) a n d c o n s i d e r t h e two ex- tremely simple situations.

1 ) t h e p e s t s f e e d only on t h e "young" trees (undergrowth);

2 ) t h e p e s t s f e e d only o n t h e "old" (adult) t r e e s .

Assume t h a t in t h e a b s e n c e of food t h e p e s t density d e c l i n e s exponentially and t h a t f o r e s t - i n s e c t i n t e r a c t i o n s c a n b e d e s c r i b e d by b i l i n e a r t e r m s as in t h e case of p r e d a t o r - p r e y system models (e.g., May, 1981; Bazykin, 1985).

Thus, f o r t h e case w h e r e t h e p e s t f e e d s o n undergrowth w e o b t a i n t h e follow- ing equations:

while f o r t h e case w h e r e t h e p e s t f e e d s o n a d u l t trees

H e r e z i s i n s e c t d e n s i t y , r i s t h e mortality r a t e of t h e i n s e c t , a n d t h e t e r m s with z z and yz r e p r e s e n t t h e i n s e c t - f o r e s t i n t e r a c t i o n .

T h e g o a l of t h i s p a p e r i s t h e c o m p a r a t i v e a n a l y s i s of models (A.O), (A.l) a n d (A.2). In t h e final p a r t of t h e p a p e r w e c o n s i d e r biological implications of t h e r e s u l t s a n d outline possible d i r e c t i o n s f o r e l a b o r a t i n g t h e model. The main t o o l s f o r our investigation are t h e b i f u r c a t i o n t h e o r y of dynamical systems a n d t h e nu- m e r i c a l methods of t h i s t h e o r y .

2. R d t s of the investigation of model (kl)

By a l i n e a r c h a n g e of v a r i a b l e s , p a m m e t e r s , a n d time, t h e system (A.l) c a n b e t r a n s f o r m e d i n t o t h e form:

where t h e p r e v i o u s notations are p r e s e r v e d f o r new v a r i a b l e s and p a r a m e t e r s which h a v e t h e same s e n s e as in system (0.1). The new p a r a m e t e r s c a n b e p r e s e n t e d in terms of t h e old o n e s as:

In t h e f i r s t o c t a n t (i.e. w h e r e t h e v a r i a b l e s t a k e on biologically possible values)

system (1.1) c a n h a v e from o n e t o f o u r equilibria. The o r i g i n , E o

=

( 0 , 0 , 0 ) , is always a n equilibrium point. On t h e invariant plane z

=

0 , where t h e system coin- c i d e s with system (0.1), e i t h e r one o r two equilibria with nonzero c o o r d i n a t e s may e x i s t . A s in system (0.1), t h e two equilibria E l

=

( z l , y l , O ) and

E 2

= ( z 2 , y 2 , 0 ) where

a p p e a r in system (1.1) on t h e line:

On t h e line

equilibrium E l c o a l e s c e s with equilibrium E,, and d i s a p p e a r s from

R : .

Besides t h e e q u i l i b r i a E l , j =0,1,2, system (1.1) could h a v e a n additional equilibrium

* : means that new variables were introduced but, for the sake of simplicity, the old notations were preserved :

A,,

a 2 b 4

This equilibrium a p p e a r s in

R :

f o r p a r a m e t e r v a l u e s ( ~ , h ) falling to t h e r i g h t of t h e line:

in t h e p a r a m e t r i c p o r t r a i t (Figure 3). Eg p a s s e s t h r o u g h t h e p l a n e z = 0 a n d coalesces o n t h i s plane with e i t h e r equilibrium El o r Ez (Figure 4). Line S i s t a n g e n t to line D l at point

in t h e ( p , h ) - p l a n e . Line S i s divided by point M into t w o p a r t s : S1 a n d SZ on which equilibrium E 3 collides with e i t h e r El or E 2 , r e s p e c t i v e l y .

In addition to t h e s e b i f u r c a t i o n s of t h e e q u i l i b r i a , autooscillations (i.e. neu- t r a l l y s t a b l e oscillations) c a n "emerge" a n d "vanish" in system (1.1). T h e s e e v e n t s t a k e p l a c e o n l i n e s R a n d P o n t h e p a r a m e t e r plane, while t h e autooscillations ex- i s t in r e g i o n s 5 a n d 6.

Equilibrium E 3 loses i t s stability o n line R d u e to t h e t r a n s i t i o n of t w o com- plex conjugated eigenvalues f r o m t h e l e f t to t h e r i g h t half of t h e complex plane.

This s t a b i l i t y c h a n g e r e s u l t s in t h e a p p e a r a n c e of a s t a b l e limit c y c l e in system (1.1) (Andronov-Hopf b i f u r c a t i o n ) .

T h e r e i s also a line c o r r e s p o n d i n g to d e s t r u c t i o n of t h e limit c y c l e s : line P o n t h e ( p , h ) - p l a n e . On l i n e P, a s e p a r a t r i x c y c l e formed by outgoing s e p a r a t r i c e s of s a d d l e s El a n d Ez e x i s t s (Figure 5). As t h e system a p p r o a c h e s line P in parameter s p a c e (Figure 3), t h e p e r i o d of t h e limit c y c l e i n c r e a s e s to infinity, a n d at t h e c r i t - i c a l p a r a m e t e r value, t h e limit c y c l e coalesces with t h e s e p a r a t r i x c y c l e a n d d i s a p p e a r s .

The point M p l a y s a k e y role in t h e p a r a m e t r i c plane. This point i s a common point for a l l b i f u r c a t i o n lines: S I , S 2 , D 1 , D 2 , R a n d P. I t c o r r e s p o n d s to t h e ex- i s t e n c e of a n equilibrium with t w o z e r o eigenvalues in t h e p h a s e s p a c e of t h e sys- tem. This f a c t allows u s to p r e d i c t t h e e x i s t e n c e of lines R a n d P.

F o r p a r a m e t e r v a l u e s c l o s e to t h e point

M

t h e r e i s a two-dimensional s t a b l e - c e n t e r manifold in t h e p h a s e s p a c e of system (1.1) o n which all e s s e n t i a l bifurca- t i o n s t a k e place. The c e n t e r manifold i n t e r s e c t s with i n v a r i a n t p l a n e z =O along a c u r v e . Thus w e h a v e a dynamical system o n t h e two-dimensional manifold with t h e s t r u c t u r a l l y u n s t a b l e equilibrium with t w o z e r o eigenvalues a n d t h e i n v a r i a n t c u r v e . This b i f u r c a t i o n h a s b e e n t r e a t e d in g e n e r a l form by Gavrilov (1978) in con- nection with a n o t h e r problem. I t was shown t h a t t h e only l i n e s originating in point M are t h e b i f u r c a t i o n lines mentioned above.

The locations of t h e R a n d P l i n e s were found numerically o n a n IBM-PC/XT compatible c o m p u t e r with t h e h e l p of s t a n d a r d p r o g r a m s f o r computation of c u r v e s (Balabaev a n d Lunevskaya, 1978). The additional a s s o c i a t e d numerical p r o c e d u r e s a r e d e s c r i b e d in t h e Appendix. W e a l s o used a n i n t e r a c t i v e p r o g r a m f o r t h e in- t e g r a t i o n of o r d i n a r y d i f f e r e n t i a l e q u a t i o n s

-

PHASER (Kocak, 1986). F i g u r e s 6, 7, a n d 8 show t h e c h a n g e s in system b e h a v i o r as i n c r e a s e s in h move t h e system t h r o u g h r e g i o n s 3, 6, a n d 7.

3. Results of the investigation of model ( k 2 )

Model (A.2), which r e p r e s e n t s a p e s t a t t a c k i n g exclusively old trees, c a n b e t r a n s f o r m e d by scaling i n t o t h e following form:

w h e r e t h e meaning of v a r i a b l e s a n d p a r a m e t e r s i s t h e same as in system (1.1).

System ( 2 . 1 ) c a n h a v e from o n e t o f o u r equilibrium points in t h e f i r s t o c t a n t

R :

: Eo

=

( 0 , 0 , 0 ) , E l

=

( z l , y 1 , 0 ) , E 2

=

( z 2 , y 2 , 0 ) , a n d E 3

=

( z ~ , Y ~ , z ~ ) . Equilibria E l and E 2 o n t h e i n v a r i a n t p l a n e z

=

0 h a v e t h e same c o o r d i n a t e s as in system ( 1 . 1 ) ; t h e y a l s o b i f u r c a t e in t h e same manner o n l i n e s Dl a n d D 2 . AS in system ( 1 . 1 ) , t h e r e i s a n equilibrium point of system ( 2 . 1 ) in

R :

:

This equilibrium a p p e a r s in

B :

below t h e line

But equilibrium E 3 d o e s n o t l o s e i t s stability so autooscillations in system ( 2 . 1 ) are not possible. Figure 9 shows t h e p a r a m e t r i c p o r t r a i t s of system ( 2 . 1 ) . The re- gion numbers in Figure 9 c o r r e s p o n d t o t h o s e in Figure 4.

Consider in more d e t a i l t h e system b e h a v i o r in p a r a m e t e r r e g i o n 3 w h e r e damped oscillations are possible. In t h e a b s e n c e of p e s t s (i.e. z =O ) t h e system t e n d s t o equilibria E 2 with c o n s t a n t d e n s i t i e s of "young" a n d "old" trees. If a small number of p e s t s t h e n invades t h e f o r e s t , a n o u t b r e a k o c c u r s a n d t h e system moves t o e q u i l i b r i a E 3 with lower tree densities and a low density i n s e c t population. The maximum i n s e c t density r e a c h e d during t h e o u t b r e a k e x c e e d s t h a t of equilibrium

A potentially u n e x p e c t e d system b e h a v i o u r c a n o c c u r if t h e system i s at equilibrium E 3 b u t t h e p e s t density t h e n declines, p e r h a p s d u e to p e s t c o n t r o l o p e r a t i o n s or t h e influence of random environmental variation. A new p e s t out- b r e a k r e s u l t s (Figure 1 0 ) . T h e r e f o r e , random declines in p e s t density may r e s u l t in r e p e a t e d o u t b r e a k s .

4. P a r a m e t e r e s t i m a t i o n f o r m o d e l ( k 2 )

Our goal h e r e is t o d e m o n s t r a t e how t h e model (A.2) might b e applied t o a r e a l f o r e s t - p e s t ecosystem. This could lead to insight a b o u t t h e dynamics of t h e ecosys- tem o r to a determination of t h e r a n g e of applicability of t h e model f o r d e s c r i b i n g ecosystem dynamics.

The eastern s p r u c e budworm-forest system was picked as a n a p p r o p r i a t e can- d i d a t e b e c a u s e of t h e availability of s u i t a b l e information f o r many p a r t s of t h e model, b e c a u s e of t h e similarity of t h e main model f e a t u r e s with some k e y a s p e c t s of t h e budworm-forest system, a n d because p r e v i o u s models (e.g., Jones 1979, S t e d i n g e r 1984) of t h e budworm-forest system have emphasized d i f f e r e n t elements (e.g., foliage, i n s e c t p r e d a t o r s , i n s e c t d i s p e r s a l ) of t h i s system.

The e a s t e r n s p r u c e budworm, Choristoneura fimiJerana (Clem.), i s a n a t u r - ally o c c u r r i n g d e f o l i a t o r of balsam f i r (Abies balsamea [L.] Mill.) in t h e b o r e a l f o r e s t s of e a s t e r n North America. Outbreaking populations kill t h e i r h o s t trees o v e r wide areas. Outbreak c y c l e s r a n g e f r o m 26-40 y e a r s in length with o u t b r e a k s lasting f o r 6-15 y e a r s . During o u t b r e a k s , i n s e c t numbers c a n i n c r e a s e o v e r f o u r o r d e r s of magnitude in s t a n d s of m a t u r e a n d o v e r m a t u r e balsam f i r which are p a r - t i c u l a r l y v u l n e r a b l e t o a t t a c k .

In a c c o r d a n c e with t h e simplistic n a t u r e of t h e model, which r e d u c e s t h e com- p l e x budworm-forest ecosystem to a system of t h r e e d i f f e r e n t i a l equations, w e t a k e a ' b r o a d b r u s h " a p p r o a c h to p a r a m e t e r estimation. F i r s t w e identify r e a l i s t i c r a n g e s f o r t h e p a r a m e t e r values a n d t h e n w e s e l e c t from t h e r a n g e to see how well t h e model c a n simulate t h e b e h a v i o u r of t h e ecosystem.

W e begin by estimating h , t h e n a t u r a l mortality rate of old trees in equation (A2). MacLean (1985) gives t h e "annual n e t probability of n a t u r a l mortality ( b e f o r e o u t b r e a k ) " as 1-3.8% f o r balsam f i r . Hence, if na i s t h e number of trees in a c o h o r t of old trees of a g e a , t h e n na +l

=

na e ^-h , a n d . O 1 S -a +I

4 . 0 3 8

.

Hence .O1 S h S .04 y r

.

(3.3) The p a r a m e t e r

1

r e p r e s e n t s t h e aging of trees in t h e model. However, depending on how one defines "old" trees,

1

can t a k e on d i f f e r e n t values. F o r in- s t a n c e , Bakuzis and Hansen (1965) r e p o r t t h a t balsam f i r r e a c h e s s e x u a l maturity at 30-35 y e a r s ; becomes moderately susceptible to a t t a c k a t o v e r 40 y e a r s of age, and becomes very susceptible at o v e r 6 0 y e a r s of age. Moreover, s t a n d s are gen- e r a l l y 40-60 y e a r s of a g e when established seedlings f i r s t a p p e a r . Thus w e assume t h a t trees spend a mean duration of 30-70 y e a r s in t h e physiologically young a g e group. If t h i s duration h a s a n exponential distribution with a mean of 30-70 y e a r s , t h e n

The function y ( y ) d e s c r i b e s t h e dependence of t h e n a t u r a l mortality of young trees on y , t h e density of old trees. MacLean (1985) suggests t h a t n a t u r a l tree mortality might fall in t h e r a n g e . O 1

-

.04 p e r y e a r . Hence, since c

=

minimum of y ( y ), w e approximate

c

=

.OI y r - I . (3.5)

The increased mortality at low y (old tree density) could b e a s c r i b e d to competi- tion with f e n s , s h r u b s , and hardwoods (Bakuzis and Hansen 1965) invading s i t e s opened up by t h e removal of t h e f i r overstory. Competition with o l d e r trees ac- counts for t h e i n c r e a s e in young t r e e mortality at l a r g e y . Assuming t h a t t h e in- t e r s p e c i f i c competition i s much l e s s detrimental than t h e suppression by t h e o l d e r a g e group, t h e n b

<<

y,,,

.

Taking y

,,,

^W 2.471 ( i n units

01

1 0 t r e e s 3 / h a ) ,

indicating a f a i r l y good s i t e (Bakuzis a n d Hansen, 1 9 6 5 , Table go), we a r b i t r a r i l y set

b -N .I X ymax RJ .2471 ( i n units of 1 0 t r e e s 3 / ha )

.

(3.6) Then, s i n c e 7(yma,) ^a .04 (MacLean 1985),

Ymax

=

a ( y m a x + ) 2 + c a . 0 4 S u b s t i t u t i n g with (3.5), (3.6), a n d t h e n solving f o r a ,

a

=

.00606 ( i n units of ha 2 ( ~ ~ 3 t r e e s ) - 2 y r -I). (3.7) W e h a v e now estimated a l l t h e p a r a m e t e r s of t h e f o r e s t s e c t i o n of t h e model (A.2) e x c e p t p , t h e rate of production of seedlings. This p a r a m e t e r combines f e r - tility, germination rate, a n d s u r v i v o r s h i p well p a s t t h e f i r s t y e a r of life ( i . e . , i n t o t h e middle of t h e r a n g e of a g e s of t h e 'young' a g e g r o u p ) . Hence, i t i s a difficult p a r a m e t e r t o estimate.

Our a p p r o a c h i s t o solve t h e system (A.0) f o r p using r e a s o n a b l e z a n d y values f o r t h e equilibrium without p e s t s . F o r instance, y

=

0 in system (1) with z =O when z

=

y h / f . From y,,, = 2 . 4 7 1 , a n d from Bakuzis a n d Hansen (1965, Table g o ) , t h e c o r r e s p o n d i n g value of z l i e s in t h e r a n g e 4.94

-

^7.42

l o 3

t r e e s / h a . Hence, if w e c h o o s e f

=

.017 y r s a y ( a f t e r equation (4)) a n d h

=

.04 y r a f t e r (3), t h e n t h e value of z at t h e u p p e r equilibrium ( E 2 in F i g u r e 2 ) i s a p p r o x i m a t e l y

zmax

=

5.81

lo3

t r e e s / ha.

S i n c e t h i s i s a r e a s o n a b l e v a l u e of

z , , ,

(Bakuzis a n d Hansen 1 9 6 5 Table 9 0 ) w e a d o p t

f

=

.017 y r - I a n d

h

=

.04 y r as r e a s o n a b l e initial g u e s s e s for t h e s e p a r a m e t e r s .

F o r a f o r e s t equilibrium to o c c u r n e a r (z,,,, y ),,, r (5.81,2.47) r e q u i r e s t h a t t h e f i r s t equation in system (A.2) with z =O a l s o meet equilibrium conditions at t h i s point. T h e r e f o r e , using (3.5), (3.6), (3.7) a n d (3.8),

This completes t h e estimation of p a r a m e t e r s f o r t h e f o r e s t s e c t i o n of t h e model a n d l e a v e s p a r a m e t e r s c, A a n d B to b e estimated. These t h r e e p a r a m e t e r s r e p r e s e n t t h e n a t u r a l p e s t mortality a n d t h e i n t e r a c t i o n between t h e f o r e s t a n d t h e p e s t .

F i r s t w e e s t i m a t e r , t h e instantaneous rate of p e s t mortality. A f t e r a n out- b r e a k t h e r e are o f t e n few m a t u r e a n d o v e r m a t u r e balsam f i r trees l e f t . Hence, w e assume y i s small a f t e r a n o u t b r e a k , so t h e p e s t equation in model (A.2) becomes approximately 2 N -&z. This equation h a s t h e solution zt zt r e

-".

Thus, a f t e r comparison i t c a n b e s e e n t h a t E c o r r e s p o n d s to t h e n e g a t i v e p a r t of t h e v e r t i c a l a x i s of Royama's (1984) F i g u r e 8. From t h e minimum of h i s smooth eye-drawn c u r v e w e estimate

1 S r S 1 . 5 y r - l . (3.11)

Next c o n s i d e r A, t h e i n s t a n t a n e o u s rate of tree mortality c a u s e d p e r p e s t . During o u t b r e a k s annual budworm-caused tree mortality p e a k s at 8

-

^{1 5 X} p e r y e a r (MacLean 1985). Hence, c o n s i d e r i n g budworm-caused tree mortality in isolation,

6 =

-Ayz

.

Then, assuming

z

i s r e l a t i v e l y c o n s t a n t d u r i n g t h e p e a k of an o u t b r e a k (Royama, 1984, F i g u r e I ) , yt yt # e Hence, in analogy with t h e d e r i v a t i o n of (3.3), .08 r 1

-

^{e -A' S} .15. S i n c e z p e a k s o n t h e o r d e r of

z,,, m 2 0

lo3

l a r v a e / tree (Miller 1975), t h i s r e l a t i o n s h i p becomes:

.00417 S A S .0081 i n trees l a r v a e y r

The p e r c a p i t a rate of p e s t i n c r e a s e p e r tree, B, remains to b e estimated.

When z i s small a n d y i s n e a r i t s equilibrium density, y i s r e l a t i v e l y c o n s t a n t so t h e p e s t equation in system (A.2) gives zt + l / z t ^@ e - W - = ) . I n analogy with t h e d e r i v a t i o n of (3.11). we n o t e t h a t (By -c) c o r r e s p o n d s to t h e positive v e r t i c a l a x i s of Royama's (1984) F i g u r e 8. From t h e maximum of h i s c u r v e we e s t i m a t e

S i n c e y yma, r 2.47 and s u b s t i t u t i n g (3.11)

0.8

s

B

s

1.42 ha tree yr-I

W e t h u s a r r i v e at t h e following t a b l e of p a r a m e t e r s f o r t h e model:

Table 1.

I

^{parameter unite 'range}

^I

^{initial guess}

^{I I}

h a 2 (lo3 t r e e s ) - 2 yr-l

l o 3

^trees/ha

yr-l yr-l y r -1 y r -1 y r - I

10 -3 t r e e s l a r v a e yr 10 -3 ha tree -1 y r -1

Table 2.

initial conditions:

s t a t e variable

2 (young t r e e s ) y (old t r e e s ) Z (insect larvae)

unite value

lo3

^trees/ha

lo3

larvae/tree

5.81 2.47 .005

The r e s u l t s of model (A.2), numerically integrated by a computer, a r e presented in Figure 11. The p a r a m e t e r s and initial conditions a r e chosen in a c c o r - dance with Tables 1 and 2. I t c a n b e s e e n t h a t t h e chosen p a r a m e t e r values belong t o region 3 on t h e r i g h t p a r a m e t e r p o r t r a i t in Figure 9, s o an o u t b r e a k i s expected. Computer simulation shows t h e o u t b r e a k h a s c h a r a c t e r i s t i c s resembling a s p e c t s of r e a l f o r e s t data. The o u t b r e a k length i s about 15 y e a r s which coincides well with observations (Royama, 1984). S o t h e model, despite i t s e x t r e m e simpli- city, could r e p r o d u c e limited time s e r i e s of a r e a l o u t b r e a k and c a n b e considered a s a compressed r e p r e s e n t a t i o n of some a s p e c t s of available f o r e s t d a t a .

T h e r e are two obvious differences between t h e computed o u t b r e a k s h a p e and r e a l f o r e s t o u t b r e a k s . F i r s t , t h e time of intensive t r e e mortality i s different. In t h e model t h i s t a k e s p l a c e a t t h e peak of t h e o u t b r e a k , while in t h e f o r e s t t h e mor- tality of t r e e s comes a f t e r t h e insect peak. I t may b e t h e r e s u l t of excluding con- sideration of foliage in t h e equations. In r e a l i t y , t h e insects f i r s t defoliate t r e e s and only then d o t r e e s begin t o die due t o defoliation. Nonetheless, t h i s distinction is essentially a minor detail given t h e "broad b r u s h " treatment of t h e problem employed h e r e .

A more important problem with t h e model's behaviour as far as r e p r e s e n t i n g budworm-forest dynamics i s t h e inability of t h e modelled stand t o fully r e c o v e r a f t e r t h e initial o u t b r e a k . For instance, in simulated y e a r s 50-60, t h e density of old trees (y) p e a k s a t about 3 / 4 of i t s original (t=O) value. This behaviour (damped oscillation) i s determined by t h e model's s t r u c t u r e and p a r a m e t e r values which place t h e system (A.2) in phase p o r t r a i t 3 of Figure 4. An obvious question i s whether random variation within t h e given r a n g e s of p a r a m e t e r values (Table I ) , as might o c c u r with changes in weather f r o m y e a r t o y e a r , could occasionally move t h e system into d i f f e r e n t phase p o r t r a i t s and thus maintain t h e oscillations.

Maintenance of t h e oscillations ( p e r h a p s a s a limit cycle) might also be accom- plished by a more a c c u r a t e r e p r e s e n t a t i o n of t h e ecological p r o c e s s e s considered in model (A.2). An obvious s t a r t i n g point h e r e would b e with t h e t e r m p p . This t e r m r e p r e s e n t s t h e rate of seedling establishment as a linear function of mature t r e e density. In f a c t , although a dense o v e r s t o r y of mature trees may produce many seeds, i t c a n inhibit seedling establishment by limiting t h e available light. Hence, f o r e s t reproductivity, p, may b e b e t t e r d e s c r i b e d by a s a t u r a t i n g function of mature tree density:

Thus p ( y ) . y -+ py when y is s m a l l .

-4 z,,~ when y i s Large.

Here z,,,, a constant, i s t h e u p p e r limit t o seedling establishment when y is large.

5. Discusaion of the results

The basic model (0.1) with two a g e c l a s s e s d e s c r i b e s e i t h e r a f o r e s t approach- ing a n equilibrium state with a constant r a t i o of "young" and "old" t r e e s ( z

=

h y ), o r degradation of t h e ecosystem (and, presumably, replacement by o t h e r species).

Models (1.1) and (2.1) have regions on t h e p a r a m e t e r plane ( 0 , l and 2) in which t h e i r behavior i s completely analogous t o t h e behavior of system (0.1). In t h e s e regions t h e system e i t h e r d e g e n e r a t e s or tends to t h e s t a t i o n a r y state with z e r o p e s t density. In t h i s case t h e p e s t i s "poorly adapted" to t h e tree s p e c i e s and can not s u r v i v e in t h e ecosystem.

In systems (1.1) and (2.1) t h e r e are also regions ( 4 and 3) where t h e station- a r y f o r e s t state with z e r o p e s t density exists, but is not s t a b l e t o small p e s t "inva- sions". A f t e r a small invasion of pests, t h e ecosystem a p p r o a c h e s a new s t a t i o n a r y

s t a t e with nonzero p e s t density. The p e s t s u r v i v e s in t h e f o r e s t ecosystem.

The main qualitative d i f f e r e n c e in t h e b e h a v i o r of models (1.1) and (2.1) i s in t h e e x i s t e n c e of density oscillations in t h e f i r s t system b u t n o t in t h e second one.

This means t h a t a small invasion of p e s t s a d a p t e d to feeding upon young trees in a two-age c l a s s system could c a u s e p e r i o d i c a l oscillations in t h e f o r e s t a g e s t r u c t u r e and r e p e a t e d o u t b r e a k s in t h e number of p e s t s (i.e., z,y

,z

/ y a n d z become p e r i o d i c functions of time). I t should b e mentioned t h a t t h e e x i s t e n c e of s u c h oscil- lations i s usual f o r simple models of even-aged p r e d a t o r - p r e y i n t e r a c t i o n s .

In o u r c a s e , however, t h e "prey" i s divided into i n t e r a c t i n g a g e c l a s s e s a n d t h e " p r e d a t o r " f e e d s only o n o n e of them. I t i s t h e p e s t invasions which induce t h e oscillations in t h e r a t i o ,

z

/ y , of t h e a g e c l a s s densities. Moreover, in t h e case of model (2.1), t h e p e s t invasion c a n include damping oscillations in t h e a g e s t r u c - t u r e .

When w e move o n t h e p a r a m e t e r plane towards s e p a r a t r i x c y c l e line P, t h e amplitude of t h e oscillations i n c r e a s e s and t h e i r p e r i o d t e n d s t o infinity. The oscillations develop a s t r o n g r e l a x a t i o n c h a r a c t e r with i n t e r v a l s of slow a n d r a p i d v a r i a b l e change. F o r example, in t h e dynamics of t h e p e s t density z ( t ) t h e r e a p p e a r p e r i o d i c long i n t e r v a l s of almost z e r o density followed by r a p i d density o u t b r e a k s . Line P i s a boundary of oscillation e x i s t e n c e a n d a b o r d e r a b o v e which a small invasion of p e s t s l e a d s to complete d e g r a d a t i o n of t h e system. In r e g i o n s 7 a n d 8 a small addition of i n s e c t s to a f o r e s t system, which was in equilibrium without p e s t s , r e s u l t s in a p e s t o u t b r e a k a n d t h e n tree and p e s t extinction.

I t c a n b e s e e n that t h e introduction of p e s t s feeding only upon t h e "young"

trees dramatically r e d u c e s t h e r e g i o n of s t a b l e ecosystem e x i s t e n c e . The e x i s t e n c e becomes impossible in r e g i o n s 7 a n d 8.

W e h a v e c o n s i d e r e d t h e main dynamical regimes possible in models (1.1) a n d (2.1). B e f o r e p r o c e e d i n g , however, l e t u s discuss a v e r y important t o p i c of time

scales of t h e p r o c e s s e s under investigation. I t is well known t h a t insect p e s t dynamics r e f l e c t a much more r a p i d p r o c e s s than t h e r e s p o n s e in tree density. I t seems t h a t t h i s d i f f e r e n c e in t h e time s c a l e s should b e modeled by introduction of a small p a r a m e t e r p < U into t h e equations f o r p e s t density in systems (1.1) and (2.1):

2

+ k .

But i t c a n b e shown t h a t t h e p a r a m e t r i c p o r t r a i t s of t h e systems are r o b u s t to t h i s modification. The r e l a t i v e positions of lines D l , D z and S as well as t h e coordinates of t h e key point M depend on r a t i o t / B which is invariant under substitutions E + C / p, B + B / p. The topology of t h e p h a s e p o r t r a i t s i s not a f f e c t e d by t h e introduction of a small p a r a m e t e r p, but in t h e v a r i a b l e dynamics i n t e r v a l s of slow and r a p i d motions a p p e a r . Recall t h a t model (1.1) had oscillations of a simi- l a r relaxation c h a r a c t e r n e a r line P of t h e s e p a r a t r i x cycle without any additional small p a r a m e t e r p. SO w e could say t h a t w e have a n "implicit small p a r a m e t e r " in system (1.1).

To demonstrate potential extensions of t h i s a p p r o a c h , w e now consider some qualitative implications t h a t atmospheric change might have on forest-pest ecosys- t e m s . A s suggested by Antonovsky and K o n u k h i n (1983), a n i n c r e a s e in t h e amount of SO2 o r o t h e r pollutants in t h e atmosphere could lead t o a d e c r e a s e of the growth rate p and a n i n c r e a s e of t h e mortality rate h . Thus, i n c r e a s e in atmospheric pol- lution could r e s u l t in a slow d r i f t along some c u r v e on t h e (p,h)-plane (Figure 12).

Suppose t h a t t h e p a r a m e t r i c condition h a s moved from position 1 t o position 2 o n t h e plane b u t remains in a region (8) where a s t a b l e equilibrium c a n e x i s t without p e s t s (Figure 4). But now, if t h e system i s exposed to p e s t invasions, both t h e forest and t h e p e s t become extinct. T h e r e f o r e , slow atmospheric changes could induce both vulnerability of f o r e s t s to pests, and f o r e s t d e a t h unexpected f r o m t h e point of view of t h e f o r e s t ' s internal p r o p e r t i e s .

6. Summary

I t i s obvious t h a t both models (A.l) and (A.2) a r e e x t r e m e l y schematic.

Nevertheless, t h e y seem to b e among t h e simplest models allowing t h e complete qualitative analysis of a system in which t h e p r e d a t o r differentially a t t a c k s v a r i - o u s a g e c l a s s e s of t h e p r e y .

The main qualitative implications from t h e p r e s e n t p a p e r c a n b e formulated in t h e following, to some e x t e n t metaphorical, form:

1. P e s t s feeding on young trees destabilize f o r e s t ecosystems more t h a n p e s t s feeding on old trees. This s u g g e s t s a possible explanation of t h e common o b s e r v a t i o n t h a t in r e a l ecosystems p e s t s more f r e q u e n t l y f e e d upon old trees t h a n on young trees. P e r h a p s systems in which t h e p e s t f e e d s on young t r e e s are l e s s s t a b l e a n d more v u l n e r a b l e t o e x t e r n a l impacts t h a n systems with t h e p e s t feeding on old trees. This may h a v e led t o t h e elimination of t h e l e s s s t a b l e systems o v e r evolutionary time.

2. An invasion of a small number of p e s t s into a n existing s t a t i o n a r y f o r e s t ecosystem could r e s u l t in intensive oscillations of t h e a g e s t r u c t u r e of t h e tree population.

3. The oscillations could be e i t h e r damping o r periodic.

4. Slow c h a n g e s of environmental p a r a m e t e r s may make t h e f o r e s t v u l n e r a b l e t o previously unimportant pests.

T h e r e are a number of possible d i r e c t i o n s f o r extending t h e model. I t seems n a t u r a l t o t a k e i n t o a c c o u n t t h e following f a c t o r s :

1 ) more t h a n t w o a g e classes f o r t h e specified t r e e s ;

2 ) c o e x i s t e n c e of more t h a n o n e t r e e s p e c i e s a f f e c t e d by t h e p e s t ;

3 ) introduction of more t h a n o n e p e s t s p e c i e s having v a r i o u s i n t e r s p e c i e s r e l a - tions;

4 ) t h e r o l e of v a r i a b l e s like foliage which a r e important for d e s c r i b i n g t h e e f f e c t of defoliation by t h e p e s t ;

5 ) f e e d b a c k r e l a t i o n s between vegetation, l a n d s c a p e a n d microclimate.

Finally, w e e x p r e s s our belief t h a t c a r e f u l a n a l y s i s of simple n o n l i n e a r ecosystem models will l e a d to a b e t t e r understanding of real ecosystem dynamics a n d to a b e t t e r a s s e s s m e n t of possible environmental impacts.

Appendix: Nl~merical procedures for the bifurcation lines R and

P

1. Andronov-Hopf bifurcation lineR

.

On t h e ( p , h ) - p l a n e t h e r e i s a b i f u r c a t i o n line R along which system (1.1) h a s a n equilibrium with a p a i r of p u r e l y imaginary eigenvalues AlV2

=

*i o (Ag

<

0). I t is convenient to c a l c u l a t e t h e c u r v e R for o t h e r fixed p a r a m e t e r v a l u e s as a p r o - jection o n t h e ( p , h ) - p l a n e of a c u r v e

r

in t h e d i r e c t p r o d u c t of t h e p a r a m e t e r plane by p h a s e s p a c e

R :

(Bazykin e t a l . , 1985). The c u r v e l7 in t h e 5-dimensional s p a c e with c o o r d i n a t e s ( p , h

,z

,y , z ) i s determined by t h e following system of alge- b r a i c equations:

w h e r e G i s a c o r r e s p o n d i n g Hurwitz d e t e r m i n a n t of t h e l i n e a r i z a t i o n m a t r i x

E a c h point o n c u r v e

r

implies t h a t at p a r a m e t e r values ( p , h ) a point ( z , y , z ) i s an equilibrium point of system (1.1) ( t h e f i r s t t h r e e e q u a t i o n s of (*) are s a t i s f i e d ) with eigenvalues AlV2

=

*i o ( t h e last equation of ( 8 ) i s satisfied).

One point on t h e c u r v e l? is known. It corresponds t o point M on t h e p a r a m e t e r plane at which system (1.1) has t h e equilibrium 1 0 with

A 1 = A 2 =

0 (e.g..

B

*io

=

0). Thus, t h e point

lies on c u r v e

r

and c a n b e used as a beginning point f o r computations. The point- by-point computation of t h e c u r v e w a s done by Newton's method with t h e help of a standard FORTRAN-program CURVE (Balabaev and Lunevskaya, 1978).

2. Separatrix cycle line P

.

Bifurcation line P on t h e p a r a m e t e r plane w a s a l s o computed with t h e help of program

CURVE

as a c u r v e where a "split" function F f o r t h e s e p a r a t r i x connect- ing saddles

E2,1

vanishes:

F ( p , h )

=

0.

For fixed p a r a m e t e r values this function can be defined following Kuznetsov (1983). Let

W;

be t h e outgoing s e p a r a t r i x of saddle

E 2

(the one-dimensional unstable manifold of equilibrium

E 2

in

R : ) .

Consider a plane 2

=

d , where d is a small positive number; note t h e second intersection of

w2+

with this plane (Figure 13). Let t h e point of intersection b e X. The two-dimensional s t a b l e manifold of sad- dle

E l

i n t e r s e c t s with plane z

=

d along a curve. The distance between t h i s c u r v e and point X , measured in t h e direction of a tangent v e c t o r to t h e unstable manifold of

E l ,

could b e taken as t h e value of

F

f o r given p a r a m e t e r values. This function i s well defined n e a r i t s z e r o value and i t s vanishing implies t h e existence of a s e p a r a - t r i x cycle formed by t h e saddle s e p a r a t r i c e s .

For numerical computations s e p a r a t r i x

w2+

w a s approximated n e a r saddle

E 2

by i t s eigenvector corresponding t o

A1 >

0. The global p a r t of

~ 2 f

w a s defined by

t h e Runge-Kutta numerical method. Point X was calculated by a l i n e a r interpola- tion. The s t a b l e two-dimensional manifold of E l was approximated n e a r s a d d l e E l by a t a n g e n t plane, a n d a n affine c o o r d i n a t e of X in t h e eigenbasis of El was t a k e n f o r t h e value of s p l i t function F.

The initial point o n t h e s e p a r a t r i x h a s z

=

0.005. The plane z

= 6

was defined by

6 =

0 . 1 a n d t h e i n t e g r a t i o n a c c u r a c y w a s p e r s t e p . The initial point on P w a s found t h r o u g h computer experiments. A family of t h e s e p a r a t r i x c y c l e s c o r r e s p o n d i n g t o points o n c u r v e P i s shown in Figure 14.

F i g u r e 15 p r e s e n t s a n a c t u a l p a r a m e t r i c p o r t r a i t of system (1.1) f o r s = B = I , & = 2 .

REFERENCES

Antonovsky, M.Ya. (1975), Impact of t h e f a c t o r s of environment on dynamics of population (Mathematical models). pp. 218-230 in Comprehensive a n a l y s i s of the Environment Proceedings of Soviet-American S y m p o s i u m , T b i l i s i , March 25-29, 2974. Hydromet, Leningrad, 1975.

Antonovsky, M.Ya. a n d M.D. Korzukhin (1983), Mathematical modelling of economic and ecological-economic p r o c e s s e s . pp. 353-358 in Integrated glob& moni- t o r i n g of e n v i r o n m e n t a l pollution. R o c . of 11 I n t e r n . S y m p . , m i l i s t ,

U S R ,

1981. Leningrad: Gidromet.

Antonovsky, M.Ya., W. Clark and Yu.A. Kuznetsov (1987), The i n f l u e n c e of pests o n forest age s t r u c t u r e d y n a m i c s : m e simplest mathematical models. WP-87-

70. I n t e r n a t i o n a l I n s t i t u t e f o r Applied Systems Analysis, Laxenburg, Austria.

Bakuzis,

E.V.

a n d H.L. Hansen. (1965), Balsam f i r . University of Minnesota, Min- neapolis, 445p.

Bazykin, A.D. a n d F.S. Berezovskaya (1979)' Allee's e f f e c t , low c r i t i c a l population density a n d dynamics of p r e d a t o r - p r e y system. pp. 161-175 in Problems of ecological m o n i t o r i n g a n d ecosystem modelling, v . 2 . Leningrad: Gidromet (in Russian).

Bazykin. A.D. (1985), Mathematical B i o p h y s i c s of Interacting P o p u l a t i o n s . Mos- cow: Nauka (in Russian).

Bazykin, A.D., Yu.A. Kuznetsov a n d A.I. Khibnik (1985)' m r c a t i o n d i a g r a m s of p l a n a r d y n a m i c a l s y s t e m s . R e s e a r c h Computing C e n t e r of t h e USSR Academy of S c i e n c e s , Puscheno, Moscow r e g i o n (in Russian).

Balabaev, N.K. a n d L.V. Lunevskaya (1978),

Computation of a c u r v e

i n n-

d i m e n s i o n a l space. F O R T ' Softurare Series, i.2.

R e s e a r c h Computing C e n t e r of t h e USSR Academy of S c i e n c e s , Pushchino, Moscow r e g i o n (in Rus- sian).

Fleming, R.A., M.Ya. Antonovsky, Yu.A. Kuznetsov (1987),

m e response of t h e bal- Sam fir forest to a s p r u c e budworm i n v a s i o n :

A

simple d y n a m i c a l model.

WP-87-71. I n t e r n a t i o n a l I n s t i t u t e f o r Applied Systems Analysis, Laxenburg, Austria.

Jones, D.D. (1979), The b u d w o m s i t e model. pp. 91-150 in G.A. Norton a n d C.S.

Holling (eds.),

Pest Management.

Pergamon, Oxford. 350p.

Gavrilov, N.K. (1978), On b i f u r c a t i o n s of a n equilibrium with o n e z e r o a n d p a i r of p u r e imaginary eigenvalues. pp. 33-40 in

Methods of q u a l i t a t i v e t h e o t y of d i f l e r e n t i a l e q u a t i o n s .

Gorkii: S t a t e University (in Russian).

Kocak, H. (1986),

m e r e n t i a l a n d d w e r e n c e e q u a t i o n s t h r o u g h computer e t p e r i m e n t s .

New York: Springer-Verlag.

Korzukhin, M.D. (1980), Age s t r u c t u r e dynamics of high edification ability tree population. pp. 162-178 in

Problems of ecological m o n i t o r i n g a n d ecosys-

tem modelling,

3. (in Russian).

Kuznetsov, Yu.A. (1983),

One-dimensional i n v a r i a n t m a n u o l d s of ODE-systems depending u p o n parameters. FORTRAN Software Series, i.8.

R e s e a r c h Com- puting C e n t e r of t h e USSR Academy of S c i e n c e s (in Russian).

MacLean, D.A. (1985), E f f e c t s of s p r u c e budworm o u t b r e a k s o n f o r e s t growth a n d yield. pp. 1 4 8 in C.J. S a n d e r s , R.W. S t a r k , E.J. Mullins a n d J . Murphy ( e d s . ) .

Recent Advances i n S p r u c e B u d w o r m s Research,

Can. For. S e r v . , Ottawa.

May, R . M . ( e d . ) (1981),

m e o r e t i c a l Ecology. P r i n c i p l e s a n d Applications.

2nd Edition. Oxford: Blackwell S c i e n t i f i c Publications.

Miller, C.A. (1975),

S p r u c e budworm: how i t l i v e s a n d w h a t i t does.

F o r . Chron.

51:136-138.

Royama, T. (1984), Population dynamics of t h e s p r u c e budworm.

C h o r i s t o n e u r a pumiferana.

Ecological Monographs 54:429-462.

S t e d i n g e r , J.R. (1984),

A s p r u c e budworm-forest model a n d i t s a p p l i c a t i o n s for s u p p r e s s i o n programs.

F o r e s t S c i . 30:597-615.

ACKNOWLEDGEMENT

The a u t h o r s are g r a t e f u l to M. Weinreich f o r h e r a s s i s t a n c e in t h e p r e p a r a - tion of t h i s p a p e r .

Figure 1. The dependence of "young" t r e e mortality on the density of "old" t r e e s .

Figure 2. The parametric portrait of system (0.1) and relevant phase portraits.

Figure 3. The parametric portrait of system (1.1).

Figure 4 . The phase portraits of system (1.1).

Figure 5. The separatrix c y c l e in system (1.1).

Figure 6. The behavior of system (1.1): s

=

b

=

1, r

=

^{2. p}

=

6 , h

=

2 (region 3). The Y-axis extends vertically upward from t h e paper.

Figure 7. The behavior of system (1.1): s

=

b

=

1, c

=

2, p

=

6 , h

=

3 (region 6 )

Figure 8. The behavior of system (1.1): s

=

b

=

1, E

=

2. p

=

6 , h

=

3.5 (region

3.

Figure 9. The parametric portraits of system (2.1).

Figure 10. A small d e c r e a s e in t h e p e s t density may result in an insect population outbreak.

1 ^{8 . ~ 4} ,

1 . 1 18.1 Y . l 1 . 1 ..I

I Irnj

Figure 11. An outbreak time equation.

Figure 12. The p r o b a b l e parameter d r i f t under SO2 i n c r e a s e .

Figure 13. The s e p a r a t r i x s p l i t function.

Figure 14. The separatrix c y c l e s in system (1.1).

-

B I F U R C R T I O N C U R V E S 1 S = B = l E = 2

Figure 15. A commuted parametric portrait of system (1.1).